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Abstract 



This is the third in a series of papers aimed at developing a practical time-domain method for self-force 
calculations in Kerr spacetime. The key elements of the method are (i) removal of a singular part of the 
perturbation field with a suitable analytic "puncture" , (ii) decomposition of the perturbation equations in 
azimuthal (m-)modes, taking advantage of the axial symmetry of the Kerr background, (iii) numerical evo- 
lution of the individual m-modes in 2+1-dimcnsions with a finite difference scheme, and (iv) reconstruction 
of the local self-force from the mode sum. Here we report a first implementation of the method to compute 
the gravitational self-force. We work in the Lorenz gauge, solving directly for the metric perturbation in 
2+1-dimensions. The modes m = 0, 1 contain nonradiative pieces, whose time-domain evolution is hampered 
by certain gauge instabilities. We study this problem in detail and propose ways around it. In the current 
work we use the Schwarzschild geometry as a platform for development; in a forthcoming paper — the fourth 
in the series — we apply our method to the gravitational self-force in Kerr geometry. 
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I. INTRODUCTION 



The study of self-forces in curved spacetimes is enjoying a surge of activity, motivated in part 
by the tantalising prospect of detecting gravitational-wave radiation generated by compact bodies 
in orbit around black holes [T] [2]. The roots of the self- force program can be traced back to 
Dirac's treatment of radiation reaction in classical electromagnetism [3], and its extension to curved 
spacetimes by DeWitt and Brehme [2). Contemporary interest was ignited in 1997 with the first 
derivations of an expression for the gravitational self- force (GSF), now known as the MiSaTaQuWa 
formula 016]. In subsequent years, there has been careful work on the foundations of the theory 
[THE], and on various practical schemes for numerical computations |16H36| (see reviews |37 | 138] ). 

The self-force program is broadly motivated by the two-body problem in General Relativity, and, 
more specifically, by the Extreme Mass Ratio Inspiral (EMRI) system in astrophysics [2, 39j HQ] , in 
which a compact object of mass [i (e.g., a stellar-mass black hole or neutron star) is in a strong-field 
bound orbit around a massive black hole of mass M, such that the mass ratio r\ = [ijM is very 
small. This scenario is difficult to handle in Numerical Relativity (but see, e.g., |41H44j for recent 
progress). On the other hand, it is tailor-made for black hole perturbation theory |45H49| . which 
is based around expansions in powers of the mass ratio r]. The compact mass generates a metric 
perturbation (MP) h a g at order 0(rj), which in turn leads to a back-reaction on the motion of 
the compact mass. This may be interpreted as a self-force F® l{ {cx 0(r] 2 )], which acts to deflect 
the motion away from a geodesic of the background black hole spacetime. Two key issues arise 
naturally. The first is that of regularization: how may a meaningful GSF be obtained from a metric 
perturbation which diverges at the location of the compact mass? The second is that of gauge: 
which predictions of the theory are independent of the choice of gauge of the MP? Both issues have 
been addressed in some depth by the community (see Refs. [5j HH [50H55] and [561461] ) . and 
the self-force program has undoubtedly now come of age. 

A significant milestone was passed in 2008, with the first comparison of conservative GSF 
effects with Post-Newtonian (PN) theory [56] , and the comparison between two independent GSF 
implementations, based on distinct gauges [57]. Comparison of strong-field GSF effects with Post- 
Newtonian (PN) theory |62H65] and with Numerical Relativity (NR) |66] have also been made. 
Some recent highlights in the self-force program include (for instance) using the GSF for calibrating 
the Effective One-Body formalism |65[ [6?ll69| and other syntheses |70H72] ; the first long-term 
evolutions of self- forced orbits |32[ 173] ; formulations of the second-order-in-mass-ratio problem [53f - 
[55] [T4Tl76| ; and the study of the implications of self- force for the cosmic censorship hypothesis 
|77H80] . Nevertheless, and despite rapid progress, all calculations of the GSF in the literature are 
rooted to the non-rotating (Schwarzschild) black hole background. A practical formulation for 
computing the GSF on the rotating black hole spacetime, represented by Kerr's solution, is now a 
high priority for the community. 

This paper is the third in a series aimed at developing a practical time-domain scheme for Kerr 
GSF calculations. Our approach is based on m-mode regularization in the time domain, a method 
first proposed in Refs. [20] [8Tj . In Papers I [82] and II [83] of this series, we established a proof-of- 
principle by applying the method to a simpler toy model, the problem of the scalar-field self-force 
on Schwarzschild and Kerr spacetimes. In the present work we apply the method to compute the 
GSF itself, and we take the opportunity to validate our results against well-known Schwarzschild 
results. Details of a GSF calculation on Kerr will follow in Paper IV. 

Our method combines several steps: (i) formulation of the linearized Einstein equations as a Z4 
system with constraint damping, using a (generalized) Lorenz gauge, (ii) removal of a certain singu- 
lar part of the MP with a suitable puncture; (iii) the decomposition of the system into azimuthal 
(m-)modes; (iv) numerical evolution of the resulting 2+1D equations [one set of (generally) 10 
coupled equations for each m-mode] with a finite difference scheme, and finally (v) reconstruction 
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of the GSF from convergent sums over m-modes. 

In line with the original philosophy, the aim is to develop a scheme in which the regularization 
of the MP is conducted in the Lorenz gauge. A key advantage is that, in this gauge, the singu- 
lar part of the MP, which informs our choice of puncture, is (in some sense) "isotropic" in the 
vicinity of the compact mass, and consequently the regularization procedure is well-understood. 
On the Schwarzschild spacetime, the linearized Einstein equations for the Lorenz-gauge MP are 
fully separable into Fourier-harmonic modes (using tensor spherical harmonics) [35] • As far as we 
know this is not possible in Kerr spacetime, our ultimate goal. Hence a key motivation for working 
in the frequency domain has been lost: the Lorenz-gauge system cannot be reduced to a set of 
ordinary differential equations. Time-domain schemes offer distinct advantages (as explained in 
Papers I and II), as they are well-suited to the study of highly-eccentric or unbound orbits, and 
they allow the orbital evolution of self-forced orbits to proceed in a self-consistent way |73j (unlike 
frequency-domain schemes, which must resort to use a "quasi-equilibrium" approach, based on 
osculating geodesies [321 184) ). Furthermore, the linearized Einstein equations in Lorenz gauge pro- 



vide a natural starting point for any time-domain scheme, because they form a Z4 system (Sec. II C 

[ESI ESJ 



and Ref. [85, 86J), which is hyperbolic in character. 

Notwithstanding, a frequency-domain approach to the problem is also being developed, with 
significant progress over past few years. In the approach of Shah et al. \30\ 131] . the MP in a 
(modified) radiation gauge is reconstructed from Hertz potentials obeying ordinary differential 
equations, which offers significant computational advantage. This approach has recently led to a 
first computation of a GSF effect in Kerr spacetime |36j. It is hoped that the same method could 
eventually allow computation of the GSF itself, for generic bound orbits. 

The Kerr spacetime is axisymmetric, and hence it is natural to decompose the MP into azimuthal 
m-modes. This decomposition has two key advantages: firstly, it reduces the computational burden 
for the time domain scheme, which now proceeds in 2+1 dimensions (2+1D) rather than in 3+1D. 
Secondly, it allows us to consider the m = and m = 1 modes separately from the rest of the system. 
These modes contain non-radiative degrees of freedom that require more careful consideration of 
initial data and conservation laws. In addition, in our Lorenz gauge formulation it turns out that 
the m = and m = 1 modes are susceptible to gauge instabilities which disrupt the numerical 
evolutions. Controlling these gauge modes is a theme of this work, and it motivates the use of 
a generalized version of the Lorenz gauge. The cost of m-mode decomposition is two-fold: the 
m-mode versions of the puncture (and effective source) are found by evaluating integrals, which 
can be computationally costly, and the self- force must be reconstructed from a sum over modes. 
The latter does not seem to pose a practical problem, as the mode sum converges rather rapidly 
with currently available punctures (the convergence properties of the m-mode sum were carefully 
investigated in Ref. [20] and in Papers I and II). 

The remainder of this paper is organized as follows. In Sec. [H] we present a preview of the 
fundamentals of our approach, giving details of the linearized Einstein equations (II B), the Z4 
scheme in (generalized) Lorenz gauge with constraint damping ( II C ) , and the m-mode regulariza- 
tion scheme (II D). This section also includes a discussion of relevant conservation laws (Sec. HE), 



which, though presented elsewhere [87] , may be somewhat unfamiliar to many in the self- force 
community. In Sec. Ill we describe some features of our implementation for circular orbits on 
Schwarzschild spacetime. Here we give explicitly the m-mode field equations (IIIB), the m-mode 
decomposition of the puncture and effective source ( III C ) , the physical boundary conditions (HID), 
and the details of the numerical implementation (HIE). In Sec. IV we present a selection of ini- 
tial numerical results (IV A), validate the results of our code for the modes m > 2 (IV B), and 
describe the manifestation of gauge mode instabilities in the m = and m = 1 modes (IV C). 
In Sec. [V] we study the low multipoles analytically, and present closed-form solutions for certain 
Lorenz-gauge modes that we diagnose as implicated in our numerical instabilities. In Sec. VI we 
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describe two practical methods for mitigating these gauge-mode instabilities. The first involves 
the use of a particular generalized Lorenz gauge, which stabilizes the m = system, and which 
reduces to the Lorenz gauge at late time. The second is a frequency filter for the m = 1 mode, 
which eliminates the undesirable gauge modes. In Sec. VII we present results for the total GSF, 
and compare with the literature. We conclude in Sec. VIII with a discussion of the way forward 
towards Kerr calculations. The appendices contain some technical details of our calculations. 
Throughout this work we use geometrized units, with G = c = 1. 



II. FUNDAMENTALS 



In this section we describe the elements of our approach. Subsection II A is a recap of the 



bare essentials of GSF formulation. In Sec. II B we state the field equations, and in Sec. II C 
we describe the Z4 formulation with gauge constraint damping that underpins our time-domain 



approach. Subsection II D describes our puncture scheme and its m-mode implementation. Finally, 
in Sec. |II E we show that the linearized Einstein equations give rise to one quasilocal conservation 
law for each Killing vector that the (Ricci-flat) background admits. It seems that this construction, 
which can help us to understand the physical content of the non-radiative multipoles, has not been 
widely appreciated in the self- force literature so far. 



A. Gravitational self-force formulation 

The aim of the GSF program is to provide an effective description of the general-relativistic 
motion of a "small" compact body, of mass [i and characteristic size l s ~ /i, on a dominant 
background spacetime with typical length scale 1Z 3> /i. For an EMRI, expansion in the point- 
particle limit l s — > corresponds (approximately) to an expansion in powers of the small mass 
ratio r]. At leading order, the small body behaves as a point particle moving on a geodesic of 
the background spacetime. At subleading order, the point particle acts as a source for a MP 
h a /3 ~ 0(r]), which exerts its influence back on the body to generate a GSF. Some care is needed 
when handling expansions, because the MP h a $ generated by a point-particle is singular, diverging 
like ~ /x/r in the limit r — > (where r is a measure of the spatial distance to the body). Mino, 
Sasaki and Tanaka obtained the first expression for the GSF (5] by applying a method of matched 
asymptotic expansions, in which "inner" and "outer" expansions for the MP were matched in a 
buffer region defined by /i/M Cr/TJCl. The key expression was also obtained using an axiomatic 
approach by Quinn and Wald [6j. Further rigorous work has subsequently put the theory of GSF 
on a firm footing [T3] (see [3S] for an introduction to the literature). Below we state the key 
results. 

The motion of the compact body may be described by a worldline 7 on the background space- 
time. Let us parameterize this worldline by x a = z a (r), with four- velocity u a = dz a /dr, where r 
is proper time along 7. For finite rj, the motion is governed by the self-forced equation 

Ha a = Fg l ^O(r?), (1) 

where a a = u a .^uP is the self-acceleration. We use a semicolon to denote covariant differentiation 
with respect to the background spacetime, and indices are raised and lowered using the background 
metric g a p. At leading order in 77, the GSF is given by 

*3f = -f (V + u a vP) [2h%. s - hfyf) u^u 5 . (2) 
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Here h^o is the Detweiler- Whiting "R- field" [8], which is a particular locally-defined smooth vac- 
uum solution of the MP equations, obtained from the full (retarded) MP via a procedure described 



in Sec. II D below. The retarded MP itself, h a p, is a solution to the linearized Einstein equations 



(Sec. II B) with a point-particle source, i.e. 

/oo 
{-g)-y 2 5 A [x» - z v {T)]u a u p dT, (3) 
-oo 

where g is the metric determinant for the background. 

In this work we will not solve Eq. ([!]) for z a (r). Rather, we will compute the GSF along 
a fixed geodesic of the background geometry, ignoring the back reaction from the GSF on the 
motion [i.e., we assume z u (r) in Eq. ^ represents a geodesic trajectory]. We leave for future work 
the important task of implementing a self-consistent evolution scheme \73\ I84| . in which the GSF 
information is fed back at each time step in order to compute the evolving orbit. 

B. Linearized Einstein Equations 

As described above, we split the spacetime metric into a background g a p at order rp, and a 
MP h a p at order 77 , and (in the first-order formulation) we neglect higher-order terms. Further 

key quantities, such as the affine connection may be expanded in a similar way, i.e., T^ ^ + 

•qT^f +0{rj 2 ). We note in passing that, although the affine connection is not a tensor, its first-order 
variation, 

^afT = \sT ( W + Kfi;a ~ K^u) , (4) 

behaves like a tensor with respect to the background spacetime |45|, 146] . The Einstein tensor is 
expanded in a similar fashion, G afi = G% + r,G% with eg for a vacuum background. At first 
order in the mass ratio the Einstein equations take the form 

- 2G$ = A*p + B a p = -16vrT^, (5) 

where 

A af} = h a p, v v + 2R^ a 5 p h l5 , (6) 

Ba/3 = 9apZ V ;v — Z a -p — Z[S- a , (7) 

Z a = hap 13 . (8) 

Here R^"*^ a is the Riemann tensor for the background spacetime, and h a p is the trace-reversed 
MP defined by 

ha/3 = h a p — \g a ph, (9) 

with h = h^- 

A gauge transformation x a — > x a> = x a — £, a (x), where £ a ~ 0(r]), generates a MP 

fc§ = 6tf + (io) 

If the gauge vector is at least twice-differentiable, then this "pure gauge" MP is automatically a 
solution of the vacuum equations, i.e. Eq. ([5| with T a p = 0. 
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C. Gauge choice, Z4 systems and constraint damping 



The standard prescription for GSF regularization (see Sec. II D below) is formulated in the 
Lorenz gauge, defined through 

Z a = 0. (11) 

This condition leads to the simplification B^ v = 0, so that the linearized Einstein's equations ^ 
becomes hyperbolic. Within the Lorenz gauge there remains a residual gauge freedom: Consider 
a gauge vector £ a that satisfies £, a]U u = 0. Then the corresponding (trace-reversed) MP h^l = 



£a;/3 + ip-,a — 9a^ ' -u clearly satisfies the Lorenz-gauge condition ( 11 ), assuming, as we do here, that 
the background is Ricci-flat. 

In our work we shall also define a generalized version of the Lorenz gauge, 

Z a = H a (hp 1 ,x), (12) 

where H a are four "gauge-driver" functions, to be specified. In order to preserve the hyperbolic 
character of Eq. ([5]), we permit H a to depend only on the MP, and not on its derivatives. In 
addition, we insist that H a is regular everywhere, even on the worldline. The trivial choice H a = 
corresponds to the Lorenz gauge. (The idea of using gauge-driver functions is not a new one [85], 
and it is being applied in many formulations of numerical relativity.) 

A key feature of the linearized set ([5]) with Z a = is that its initial- value formulation preserves 
the Lorenz-gauge condition: If a solution has Z a = on an initial Cauchy surface, then this 
condition will be satisfied at any time. However, since suitable initial data (i.e., Lorenz-gauge data 
compatible with a moving point-particle source) are not available, and since finite-differencing 
numerical error is inevitable, in practice any numerical implementation of Eq. ^ would lead to 
gauge-condition violations. A practical solution, outlined in Refs. [13 [88], employs ideas from the 
Z4 formulation of Numerical Relativity \85\ IBS] . The Z4 system is obtained in its 4-dimensional 
covariant form by replacing the vacuum Einstein equations G a p = by G a p + V^Z^ + V ' pL a — 
0a/3 V 7 Z 7 = 0. Here Z a is a supplementary four- vector of constraints. A solution of the Z4 equations 
is a solution of the Einstein equations if and only if Z a = 0. 

In our linearized system, we take the constraint vector to be 



leading to the new system 



where A a p is as in Eq. ([5}, and 



7L a — — g [Z a — H a ) , (13) 



Aaf} + B afl = -ltorTaf,, (14) 



Ba/3 = g a pH v ' ■„ — H a -p — Hp- a . (15) 
Note that, after this replacement, the constraint vector Z Q obeys a wave equation CE a = "L a -^ = 0. 



[to see this, consider the divergence of Eq. (14) and make use of Ricci-flatness and stress-energy 
conservation T a p^ = 0.] The set (14) is only equivalent to (JH| if Z Q = (i.e., if H a = Z a ). In order 
to drive Z Q towards zero, we supplement the set (14) with a "constraint damping" term C a p |88| . 



In principle, any choice of C a p that vanishes as Z a — > is justified. We will employ a constraint 
damping term of the form 

C a/3 = -2At(n Q ,Z ( a + npL a ), (16) 
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where k and n a are scalar and vector fields to be specified in Sec. IIIB 1 

In summary, in the following sections we will implement the modified system 

■A a p + B aB + C a3 = -16irT aB , (17) 

where the left-hand terms are given by Eq. ^ , ( 15 ) and ( 16 ) . The key requirement of any successful 



numerical implementation is that the constraint violation dissipates with time, 7L a — > 0, so that 
we recover a valid solution of the linearized equations In Sees. Ill and IV, we describe a 



Lorenz-gauge (H a = 0) implementation that works well for the modes m > 2, but which suffers 
from linear-in-i gauge mode instabilities in the modes m = and m = 1. In Sec. |VI| we describe a 
generalized Lorenz-gauge implementation that goes some way towards curing the instability. 

D. Regularization 



As discussed in Sec. II A a point-like source generates a MP which diverges along the particle's 
worldline, and a method of regularization is required in order to extract the correct regular field 
h^p that enters the GSF construction formula (J2J). Detweiler and Whiting [8] gave a prescription 
for constructing h^g through a subtraction 



ha/3 - ha/3 ~ h a g, (18) 

where h aB is the retarded solution of Eq. (|5j) with the point-particle source (ph, and h^g (the S 



l a8 

h aB = u a g — lb a g, 

field, for "singular /symmetric") is locally defined in terms of a Green function G^, j3 ? l3 (x, x'), which 
is (i) symmetric in its indices and arguments, (ii) a solution of the inhomogeneous equation ([5]), 
and (iii) zero when x and x 1 are connected by a timelike geodesic. The Green function has a local 
definition in terms of a Hadamard parametrix. The terms in this parametrix can be expanded in 
powers of the coordinate separation of the points. This leads to a local expansion for the S field 

EH EH E3 ESI EDI- 

1. The puncture scheme 

Following [20] and Papers I and II, let us introduce the V field h (the "puncture"), which has 

g 

the same local expansion as the S field, h , up to a certain order, and which has a smooth global 
continuation. Corresponding to this puncture field is an 1Z ("residual") field, defined through 

= K B -h V aB . (19) 

If the puncture field agrees with the S field in the local vicinity of the worldline up to a suitably high 
order (see Paper I), then (i) the residual TZ and Detweiler- Whiting R fields agree on the worldline, 
and, (ii) the GSF can be obtained from the gradient of the TZ field using 



F^ s =^% 1 . s = (20) 

x=z(t) X-¥Z(T) 

where 

k a W(x) = \g a5 u^u^ - g^u^u 5 - \u a u^v?u s + \u a g^u s + \g aS g^. (21) 

Here u a = u a (x) is any smooth extension of the four- velocity u a off the particle's worldline. The 
so-called "red shift" variable H |31| [56| [571 191] may also be found from the residual field in a 
straightforward way: 



H EE \u a U^h% 



= hm hu a H^. (22) 

X = Z(T) X^Z(T) 
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2. m-mode decomposition 



Following the approach of |20| 181] and Papers I and II, we take advantage of the axisymmetry 
of the Kerr spacetime (i.e., the existence of the azimuthal Killing vector Xfy = dx a /d<p), to 
decompose key quantities into azimuthal m- modes. That is, we let 

oo 

X(t,r,9,<f>) = E ^ (m) (t,r,e)e im t (23) 

m=— oo 

where X is any relevant physical quantity — for example, a component of the retarded MP h a p — 
and (t, r, 6, 4>) are Boyer-Lindquist (BL) coordinates. If X is real, then it follows that X^ satisfies 
the complex-conjugation symmetry, 

X (-m) =x (m)*_ (24) 
To obtain the m- modes, we may apply the inverse transformation, 

(t, r,9) = — I" X(t, r, 6, (fye'^dcj). (25) 

2?r J-n 

Here some care is needed if X^ is not a smooth function; in particular, when applying the 
transform to the puncture field, we should first ensure that it is smooth everywhere off the particle's 
worldline. 



Key quantities, such as the value of the R-field h a a at the particle, or the GSF exerted on 
the particle, may then be recovered from a sum over modes: 

^Kr)]=EV > where Kp =e m Re lim Cj } e™>, (26) 



*5d*(r)] = £ F^\ where = ^ m Re lim W*i*V s (h% m) *°*) . (27) 

m=0 

Here e m = 2 for m ^ and e m = 1 for m = 0; we have folded over the m < contributions onto 
the m > ones, with an overhat indicating the combined contribution from the two ±m-modes 
for given m (which forms a real quantity). Some care is needed to show that reversing the order 
of the sum and limit is justified; this analysis is given in Ref. |20j. 



E. Conservation laws 

In this section we will construct quasilocal conserved quantities of the linearized system ^ 
corresponding to symmetries of the background metric. These will play an important role in our 
discussion of low multipoles in Sec. [Vj 

Let us suppose that the background spacetime admits a Killing vector X a (for the following 
discussion we only assume that the background is vacuum; we will specialize to Kerr later). Then, 
following Abbott and Deser |87| . we introduce the antisymmetric two- form 

F<*P = — g- (^X x h X [ a] /3} + X x .{ a h /3 }x + X[ a Zpjj , (28) 
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FIG. 1: Diagram illustrating the construction discussed in the text. S is a closed spacelike 3-surface, defined 
by t — const and rj < r < r 2 , outside a Kerr black hole. The hypersurface is penetrated by the particle 
worldline 7. Stokes' theorem is used to relate the particle's energy and angular momentum to surface 
integrals of the MP over the boundaries <9£i and <9£2 at r = r\ and r — r 2 , respectively. 



where Z a is defined in Eq. my and square brackets denote antisymmetrization: Yy a .^ = \{Y a -^ — 
Yg. a ). It is reasonably straightforward to verify that the divergence of F a n is proportional to the 
left-hand side of Eq. ^ contracted with the Killing vector, so that 

F a ^ = T afs X? = j a . (29) 

The divergence of the current j a is zero, j a ' a = 0, due to the conservation of stress-energy, T a/ g ;/3 = 
0, in conjunction with the Killing property X a .p + Xp- a = 0. This means that the total "charge" 

Q(X) = Jj a dZ a (30) 

contained in a spacelike hypersurface E extending to infinity is conserved (E-independent), assum- 
ing j a = at spatial infinity (see, e.g., Chap. 3 of Ref. [92J; we follow here the notation of [92J, 
with dY> a representing the appropriate vector area element on E). If the background spacetime 
is stationary, i.e., admits a time-translation Killing vector Xf t y than Q(X%\) gives a Komar-like 
definition of the mass- energy content of the MP. Similarly, if the background spacetime is axially 
symmetric, with a rotational Killing vector X^, than Q(X?U) gives a Komar-like definition of the 
angular-momentum content of the MP. [These definitions are Komar-like in that they bear on the 
time-translation and rotation symmetries of spacetime. However, here the relevant symmetries are 
these of a background spacetime; the full spacetime (background+perturbation) need not have any 
symmetry for our definitions to hold.] 

Let us now specialize to a Kerr background and a point particle with stress-energy given by 
Eq. Q. Consider a closed spacelike 3- volume E, defined by t = ts(= const) and < n < r < r%, 
for some r\ and r2 and with r = r^ being the BL radius of the event horizon (see Fig. [TJ . We find 
that the "charge" contained in E is 

Q(X) = ri<r p (t s )<r 2 , (31) 

I 0, otherwise, 

where r p (tz) it the particle's radius at time is- With the Killing vectors X^ = dx a /dt and 
X?U = dx a /d(j) of the Kerr geometry, the quantity fiX a u a on the right-hand side is, respectively, 
(minus) the energy \xE = —fiut and the angular momentum \xL = fJ,u^, which are conserved along 
the particle's geodesic worldline. 
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Using Stokes' theorem we may express Q(X) more usefully in terms of quantities on the 2- 
dimensional boundaries of E. Let dH\ and dY<2 denote the inner and outer boundaries of E, at 
r = r\ and r = r 2 , respectively. Then from Eqs. (29) and (30) it follows (cf., e.g., Sec. 3.3.3 of 
Ref. 02]) that 



Q(X) =T(X,dX 2 )-T(X,dXi), 



(32) 



where 



r{x, as) 



a/3, 



(33) 



in which dJl a is an appropriate 2-surface element on <9£ |92j . We note that the surface integral 
J-(X, dTi) (unlike F a g itself) is gauge invariant, in the sense that it vanishes for any (smooth) 

h aj3' F aP Can be 



written as the divergence, F 



(£) 



pure-gauge MP of the form (10). This can be verified by noting that, for h a t 



a/3 



n 
a/37 



, of a 3-form 



Va/3-y 



8vr 



{X[ a ^0n] + ^[a;/3? 7 ]) • 



(34) 



The integral of F^) over a 2-surface of constant t, r then works out to be proportional to f (Aq^ — 



A r f,fi)d6d(j), with As = rf-^ 'e a ^g ■ This surface integral vanishes for any A& (assuming £ c 
hence Ax, are smooth over the 2-surface), leading to J"® (X, <9£) = 0. 



and 



Eqs. (31) and (32) relate the conserved quantities £ and C to the difference between two surface 



integrals J- enclosing the relevant volume. In fact, a stronger result can be established, in the form 
of a statement about J- itself, for any 2-surface e?£ r of constant t = t% and constant r: 



F(X,d-£ r ) 



»X"u a , r>r p (t s ), 
0, r < r p (t s ). 



(35) 



Xfa or X 



X^y It follows from the following 



This result is valid in any gauge for either X 
argument. Consider the surface integral J r (X,dT loo ) for a 2-sphere 9Soo at r — > 00. Using the 
m-mode decomposition described above we may formally split the physical MP into a stationary, 
axially-symmetric piece (m = 0) and a non-axially-symmetric piece (m > 0), with the latter 
averaging to zero upon integration over <9£oo and hence not contributing to J- . The remaining, 
m = contribution is made up of two pieces, one coming from the particle's monopole-like mass 
perturbation, and the other from the particle's dipole-like angular-momentum perturbation; higher 
multipoles decay fast enough at r — > 00 for their contribution to J~(X, <9£oo) to vanish. It is 
straightforward to construct explicit expressions for the mass and angular-momentum perturbations 
in a particular gauge: Let g a a (x a ; M, J) denote the Kerr metric in BL coordinates, with M and 
J = aM being the black hole's mass and spin. Then the linear variations 



ASM) 



dM 



and 



AS J) 



dJ 



(36) 



M 



(with fixed BL coordinates) describe, respectively, the mass and angular-momentum perturbations 
associated with a particle of orbital energy \x£ and angular momentum \xC An explicit calcula- 

fi£ and -F{X^ y d^oo) = for hffi, and F^XfodEoa) = and 
Therefore, for the complete MP we find T(X, <9£oo) = ^X a u a , 



tion gives 



c 



fiC for h^p. 
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which should hold irrespective of the gauge chosen for the MP, by virtue of the gauge invariance 
of T. The result (35) then follows directly from the "jump" condition (32) with (31). 

The relation in Eq. (35) is useful in that it provides us with a simple criterion by which to 
determine if a given numerically-constructed MP has the correct mass and angular-momentum 
contents. In Sec. [V] we will show how this criterion can be used constructively to obtain the correct 
(physical) nonradiative pieces of the MP. 



III. IMPLEMENTATION: CIRCULAR ORBITS ON SCHWARZSCHILD SPACETIME 

In this section we describe the details of the first implementation of the m-mode regularization 
method for computing the GSF along geodesies of a black hole spacetime, focusing on the specific 
case of circular orbits on Schwarzschild spacetime. 



A. Setup: Circular geodesies on Schwarzschild spacetime 

In the region exterior to the black hole, the Schwarzschild geometry is described, using 
Schwarzschild coordinates {t,r,6,cj)}, by the line element 

ds 2 = g^dx^dx" = -f{r)dt 2 + / -1 (r)dr 2 + r 2 {d9 2 + sin 2 9d(f> 2 ), (37) 

where 

fir) = 1 - r h /r, (38) 

M is the mass of the black hole, and r = Th = 2M is the radius of the event horizon. To describe 
the region beyond the future event horizon % + , it is common to introduce the advanced time 
coordinate v = t + r*(r), where the "tortoise" radial coordinate r* is defined via 

r* = r + r h ln(r/r h - 1), (39) 

with dr*/dr = / _1 (r). This leads to the ("advanced") Eddington-Finkelstein line element, 

ds 2 = -f(r)dv 2 + 2dvdr + r 2 {d6 2 + sin 2 Odcf 2 ). (40) 

We consider a pointlike particle of mass \x in a circular orbit around the black hole. Neglecting 
GSF effects, the particle follows a geodesic z a (r) with a four- velocity u a = dz a /dr with respect to 
proper time r. Without loss of generality, we take the orbit to be in the equatorial plane, so that 
(in Schwarzschild coordinates) 

z a (r) = [t(r), r , tt/2, 0(t) = Qt(r)] . (41) 

The orbital frequency is 

n = d(j)/dt = Jm/t^, (42) 

and the four-velocity is 

u a = (£//o)[l,0,0,fi]. (43) 
Such a geodesic has (conserved) specific energy and angular momentum given by 

S = -Ut = f (1 - 3M/r )- 1/2 , (44) 
C = U4> = (Mr ) 1/2 (l-3M/r y 1/2 , (45) 

where /o = 1 — 2M/rQ. 
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B. 2+1D field equations in vacuum 



We next describe our method for solving the linearized Einstein equations in 2+1D, subject to 
the Lorenz-gauge constraint = (i.e. = 0). Focusing for now on the vacuum case (T a g = 0) 
and specializing to the Lorenz gauge {B a p = 0), we have the linear system of coupled PDEs, given 
by 







(46) 



[recall Eqs. ^ with Q]. 

Consider the m-mode decomposition of this system. We write, formally, 



T r( m ) imd> x( m ) M / a \ ( m ) n a\ 

Kp= 2^ Kp e > Kb = -Jap(r,8)u y a p J (t,r,9) 



(47) 



m=— oo 



(no summation over a/3), with J a p(r, 0) = g a B9aB (again, no summation) and g a p 
diag [1,/- 



, r 2 , r 2 



sin 2 0] . Here we have chosen the 7 Qj g prefactor to (i) scale the components of 
the MP in a similar way to the components of g a @, and (ii) avoid f 1 ^ 2 terms in the equations that 
follow. Inserting (47) into (46) leads to a set of ten coupled equations 

0. 



j, i — i (m) , (m) 



where = (dtu^, d r „u^ u , dgu^, u^) are given in Eqs. (A1)-(A10) of Appendix 



/□sc = 



d 2 d 2 f 

1 h — 

dt 2 dr 2 r 2 



d 2 n d 

W 2+coi0 de 



(2M m 2 
\ r 



sin 



(48) 



and 



(49) 



is a scalar- like wave operator. This set is complemented by an m-decomposed version of the four 
Lorenz-gauge equations Z a = 0, which we give explicitly in Eq. (A11)-(A14) of Appendix [A! 



1. Gauge constraint damping 



As noted above (and in Refs. [T71 [88] ) 



towards zero by adding suitable damping terms to the set (46) 
combinations of Z a in any convenient way. 
formulation to be a good choice: 



numerical gauge-constraint violations may be driven 
In principle, we can add linear 



In practice, we have found the following covariant 



(50) 



A a p + C a/3 = 0, where C Q/ g = k (n a Zp + npZ a ) , 

and k = df/dr = 2M/r 2 with n a = [1,/ -1 ,0, 0]. Note that n a is a null vector (i.e., 
and our approach is similar, but slightly different, to that taken in Refs. |17l [21~1 [58] 
that, with this choice of constraint damping, all first-order t and r* derivatives appear in Eq. (50) 
only in the combination d v = \{dt + <9 r #). With this choice, the 2+1D vacuum equations become 



n a n a = 0), 
. Note also 



t> i — i (m) , k ~a (m) 



0. 



(51) 



with M^a = ■M-^B (dvU^ujdeu^jU^) given explicitly in Eqs. (A15)-(A22) of Appendix 

Alternative choices of constraint damping are possible, such as the version used by Barack and 
Lousto [T7] (and later Barack and Sago |21j ) in the 1+1D setting. In the 2+1D setting, we find that 
implementing the Barack and Lousto choice leads to numerical evolutions which are susceptible 
to exponentially-growing instabilities at late times in modes m > 2, and at early times in modes 
m = and m = 1. 



Am) 



A 
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C. Sourced 2+1D field equations and the puncture scheme 



The stress-energy associated with a pointlike particle was given in Eq. ([3]), with (— g) 1 ^ 2 = 
r 2 sin#. After taking the m-mode decomposition (see Sec. IID2) we reach the sourced m-mode 
equations 



f^+M^ = S^, (52) 



where 



SS = -^^e- imnt 6(r - r )5(9 - Tv/2)u a up (53) 

(no summation over a/3 implied). Naively, we might attempt to solve this equation immediately 
on a 2+1D grid. However, the (5-functions are difficult to implement, and, worse, the components 
of the retarded MP, u^g , will in general diverge logarithmically as (r, 9) — > (ro,7r/2) |81j . 



Second- order puncture 



In this work, we implement a puncture scheme first introduced in Ref. [20J. There, it was shown 
that, under the m-mode decomposition, it is sufficient to use a puncture of second order (under 
the classification scheme of Paper I). This is in contrast to 3+ ID schemes, in which the third-order 
puncture is the minimum requirement to obtain a residual field which is differentiable. A heuristic 
explanation for this is given in Sec. II. G. 2 of Paper I. 

The Detweiler- Whiting S field may be defined locally for field points x in the vicinity of a fixed 
point z on the worldline. To promote the S field to a puncture function which is defined in the 
vicinity of the entire worldline, the next step is to make z a function of the field point x. The 
simplest way to do this is to choose z to be the point on the worldline with the same time t 
coordinate as x. We may then define the coordinate differences, 



Sx a = x a - z a (t) = [0,Sr,69, 



(54) 



where Sr = r — ro, 69 = 9 — ir/2 and 5(ft = (ft — fit. Now, following [20], we introduce the puncture 
function 



Kp\ x ) = — 

£ [2] 



u a Uj3 + (Tq up + Tg u a ) uxSx 1 



(55) 



where u a and are the tangent velocity and the affine connection evaluated at z a (x), and ep] 
is defined as follows: 



£ [2] 



{g a p + u a up)\ z 5x a 5x p + (uxu^T^ + g a p, 7 /2) 5x a 5x^5x 1 



(56) 



—V 

The above function h a a(x) approximates the Detweiler- Whiting S field through 0{5x ). 



There is one remaining problem: the puncture function defined in Eq. (55 ) suffers a discontinuity 



at 



±7r. To correct this, we replace the "even" and "odd" terms in 5 (ft as follows: 



b 2 -> 2(l-cos<ty) [= 5(ft 2 + 0(5(ft% 
>4> -»■ sin 5(ft [= 6(ft + 0(S(ft 3 )]. 



(57) 
(58) 
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—■p n 

This keeps the local expansion of h a g(x) near the particle intact through 0(5x ), so our puncture 
remains equal to the S field through this order. With the above replacement the puncture takes 
the form 

h ap = , (59) 

in which 

(-v = [s200^ 2 + s 2o50 2 + S3oo5r 3 + si2o5r86 2 + 2(s o2 + siojfr)(l - cos 5(/>)] 1/2 , (60) 

and the nonvanishing components of Xa/3 are 

_ f UaUfi + D a/3 5r for a/3 = tt, t<f>, <f>t, 4>4>, , s 

Xa[3 ^ D a p sin 5(f) for a/3 = tr, rt, r<f>, (f>r. 

The coefficients Sijk, which depend on tq only, were given explicitly (in the more general Kerr case) 
in Appendix A of Paper II. The coefficients D a p read 

2£ 2 M , , 

r (r - 2M) 

= D H = , (63) 

ro 
C 2 

D r(t> = Da,t = — , (65) 
2C 2 

D H = . (66) 



2. m-mode decomposition of the puncture function 

The next step is to compute the m-mode representation of the components of the puncture 



function. Using Eq. (25) with <ft = 5<f> + Qt we write 



-V(m) e T V 



h °* ~toTJ. h ap(Sr,8e,S^e-^d(Scl>). (67) 
It turns out, as in the scalar- field case explored in Papers I and II, that all the relevant integrals over 



S(j), in the puncture and source (Sec. Ill C 3 below) can be represented in terms of complete elliptic 

integrals of the first and second kinds, as we detail in Appendix |B} The m-mode representation of 
the puncture itself is straightforward: we find 

j-V{m) _ 4/i ( m ) imUt , s 

where the nonvanishing components of XaB are 

~M _ f ( u aup + Da/30~r)I™ for a/3 = tt, t(f>, <j)t, 4><j), . . 

Xa ? ~ \ DapJfi 1 for a/3 = tr, rt, r<f>, cpr, 1 > 



with I™ and J™ each being a combination of complete elliptic integrals given explicitly in Eqs. (B3) 
to (Bll) of Appendix pi Note that I™ is purely real and logarithmically divergent at the particle, 



whereas J™ is purely imaginary and it is continuous and differentiable at the particle (although 
its second derivatives diverge logarithmically there). 
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3. m-mode decomposition of the effective source 



rR 



The residual field h a g, defined in Eq. (19), obeys a wave equation with an "effective source", 



Vh 



■R 

a/3 



T 



cff 



(70) 



where here we use D as a general shorthand notation for our choice of operator on the left-hand 



side of Eq. (17), which depends on the choice of gauge and constraint damping. The effective 



source is found from the action of the operator on the puncture function, i.e., 

-v 



-7-cff 

'a/3 



-16ttT, 



a/3 



(71) 



-v 



Here, recall that is expressed as a function of {5r, 59, 5(f)}, and we may express the partial 
derivatives in T> in terms of these variables, e.g., djd<f> = d/d(5(p) and d/dt = —ft~ 1 d/d(5(j)). Let 



us now choose the left-hand side operator as in Eq. (50). The m-mode decomposed equations 
governing the residual field are then 



fa sc n^ m) +M n{m) 



"a/3 

where M^g is obtained by replacing 



a/3 
Tl(m) 



s 



eff(m) 
a/3 ' 



(72) 



from h 



K{m) 



l a/3 



in Eqs. (|Al5t)-(|A22t) , with u^ m) defined 



a/3 



as in Eq. (47). The effective source for the m-mode system is 

,eff(m) _ fr 



s : 



aj3 



7a/3 



^-eff(m) 
'a/3 



(73) 



(no sum implied), where T^™^ are the m- modes of 7^f. With our second-order puncture, the 
effective source diverges logarithmically on the particle's worldline (see Papers I & II). 

For circular orbits, the m-mode effective source may be expressed in terms of complete elliptic 
integrals, as in the scalar case. To this end it is convenient to write S a/3 
contributions, in the form 



■efifjm) 



as a sum of two 



5' 



-eff(m) 
a/3 



27T7 Q/3 



J a(3 



(m) 

a/3 



(74) 



where Z^ m ^ arises from the "scalar" part of the wave operator (i.e, the term in T> involving Dsc), 

and AZ^™^ represents all remaining terms. The first contribution can be expressed concisely in 
terms of quantities that were already computed in Paper I: we find 



J af3 



(llaUp + D a/3 5r) Y^k=l °k^k 



D a/3 Sfc=l S kJ r k 



S h IY} for a/3 = tt, tcj), <f>t, 

for a/3 = tr, rt, r<f>, 4>r, 



(75) 



where I™ and J™ are complete elliptic integrals given explicitly in Appendix [B] below, and the 
coefficients Sk (which depend on r, 9 and vq but not on m) were given explicitly in Appendix B of 
Paper I [Eqs. (B8)-(B12)]. The remaining part of the m-mode effective source may be computed 
with the assistance of a symbolic algebra package. We obtain the form 



AZ 



(m) 

a/3 



(76) 



fc=0,l,2,6,7 



fc=0,l,2 



where the two extra elliptic integrals Jgy are also given in Appendix pi [Eqs. (Bll) and (B12)], and 
the various components of the (m-independent) coefficients c k a ^ and are given in Appendix [cj 
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D. Boundary conditions 



We require that the advanced Eddington-Finkelstein (aEF) components of the MP are regular 
across the event horizon, i.e. for R = 2M and finite v, where {v = t + r*,R = r,9,(j)} are aEF 
coordinates. The aEF components are related to the Schwarzschild components via 

(77) 
(78) 

ht, (79) 

(80) 
(81) 

and h a p = h a p for all other components. We hence require that the above algebraic combina- 
tions of Schwarzschild components remain regular (in particular, finite) as r* —> — oo. The same 

combinations are also required to be regular at the level of the m-mode MP h^g(t,r, 9). 

We also require boundary conditions at the poles, 9 = 0, tt. These are determined from the 
requirement that the MP is regular (in particular, continuously differentiable) when expressed 
in coordinates which are regular at the poles (e.g., locally Cartesian coordinates based at each 

pole). This, in turn, translates to a condition on the local behavior of the m- modes h^a (t,r,9) 
for 9 — > + , 7r _ (with any fixed t, r). The form of this condition depends on the MP component in 
question. For the North pole (9 = 0) we find, for m > 0, 



y-(aEF) 
n vv 


= h a 




7>EF) 
n vR 


= htr 


- r%t 


-j-(aEF) 
n RR 




- 2f- 1 h tr + / 


r (aEF) 
n R6 


= h r e 


- r%e, 


-j-(aEF) 
n R<f> 







9 m , a/3 e {tt,tr,rr}, 

fl'"*- 1 ', a/3 S {t9,t<f>,r9,r(j)}, (82) 
0l m " 2 l, a/3 e {99,8(i 



h ( $(r,t,9)~{ 



The conditions at the South pole (9 = tt) can then be inferred from the reflection symmetry 



TQ(t,r,*-0)= -M^' <tfe{tf,r*,i rj , (83) 
I +h a p(t, r, 9), otherwise. 



E. Numerical implementation details 

1. Numerical domain and the worldtube scheme 



The puncture function ( 68 ) was constructed by truncating the series expansion of the Detweiler- 



Whiting S field at second order, and the effective source was formed by acting on the puncture 



with the relevant wave operators. Unfortunately the global behaviour of the source (74) is not 
suitable for a numerical scheme. It generally diverges in the large-r limit, and also at the poles 
(where sm.9 = 0). Vega and collaborators [331 EH E3 EH [93] avoid this problem by multiplying 
the puncture by a windowing function, to smoothly moderate the behaviour of the effective source 
far from the worldline. By contrast, and following Ref. [81j, Paper I (see Sec. IIH) and Paper II 
(see Sec. IVA), we apply here a "worldtube" scheme. 

In the 2+ ID domain, the worldtube T is defined to be a region of finite extent in r and 9 
(or order ~ M), surrounding the worldline. Outside the worldtube, we evolve the vacuum wave 
equation for the physical (i.e. retarded) MP. Inside the worldtube, we evolve the sourced equation 
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for the residual MP. Across the boundary of the tube, dl~, we convert between the two using the 
puncture field. Hence, our evolution equations consist of 

outside T, 

inside T, (84) 
across QT . 

These equations are solved on a fixed uniform grid based on t,r*,6 coordinates; see Fig. 1 in Paper 
II. Let us denote the corresponding grid spacings by At, Ar», AO. In the 2+1D domain, the 
particle's trajectory traces a straight line at = tt/2 and r* = r*o[= r*(ro)]. We lay the grid so 
that the particle's trajectory on the grid crosses through a row of grid points (of fixed r*, 0). Then 
we define a worldtube T of fixed coordinate widths {T rt , Tq} as follows: Consider a grid point with 
coordinates (t, r*,9). If (r* — r*o| < r r>fc /2 and \9 — ir/2\ < Tg/2 then the point is said to lie within 
the worldtube; otherwise it lies outside. 

In the finite difference scheme, the derivatives of the MP at a given grid point (say, point "o" ) 
are computed based on the values at several neighboring grid points. If point o lies close enough 
to dT that some of these points are "in" and others are "out" of the tube, we make the following 
adjustment. If the point o is "out", then we first demote all relevant "in" points to "out" points 
using = u^i 71 ^ + r(^7 Q ,^) _1 /i^ m ' ) , before applying the finite-difference formula. Conversely, if 
the point o is "in" , we promote all "out" points to "in" points in a similar manner. 




,,M + M^ m ) - ( 

n(rn) CjKim) 
a a(3 Jvl a/3 
(m) 



S' 



eff(m) 



a/3 



via ) i - \ _ i v 



-V(m) 



2. Finite- difference scheme: method of lines 



To evolve Eq. ( 84 ) , we used a finite difference scheme based on the well-known Method of Lines 
(see Ref. [94 1 or Sec. 4.2 of Ref. [M])- The first step is to rewrite the equations in first-order form, 
which is simply achieved by defining the auxiliary variables, 

= (85) 

Hence, for each m > we have a 40-dimensional system (from 10 real and 10 imaginary parts of 
u a p and v a p, respectively), which is first-order in time. (An exception is the m = mode, which is 
real and thus gives a 20-dimensional system.) The next step is to replace all spatial derivatives (in 
r* and 9) with their finite-difference approximants. We chose to use fourth-order finite difference 
operators, i.e. 



d 2 X 1 



dr 2 12 A r 2 
dX 1 



(-X j+2 + leXj+i - 30Xj + 16X,_i - Xj-. 2 ) , (86) 
+ 8Xj + i — 8Aj_i + Xj^2) , (87) 



<9r* 12Ar* 

where X G «^ m '}, and Xj = A(t,ro* + jAr*,9). Similar operators were used for the 9 

derivatives. 

To evolve the first-order equations (i.e., to progress forward in time t by steps of At), we applied 
a fourth-order Runge-Kutta step. For details, see Sec. IVB in Paper II. In vacuum (i.e., without 
sources), we would expect this scheme to be globally fourth-order accurate, so that doubling the 
grid resolution will reduce the finite-differencing error by a factor of ~ 2 4 . However, when the 
effective source is included, a naive application of the scheme (as here) will not be globally fourth- 
order accurate, due to the logarithmic divergence of on the worldline. 
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3. Boundary conditions 



The spatial boundaries of our grid are at r* = r*i n , r* = r* ou t and at 9 = 0, 9 = it (but see 
below how reflection symmetry allows us to replace the latter boundary with one at 6 = 7r/2). 
Taking r*- m <C — M and r* out S> M places the radial boundaries at the asymptotic domains where 
the usual "ingoing" and "outgoing" radiation conditions apply to the retarded perturbation. In 
practice, we do not actively impose these radial boundary conditions in our time-evolution scheme. 
Rather, we set the radial boundaries far enough that a sufficient portion of the late-time solution 
in the neighborhood of the particle has no causal connection with these boundaries. For the modes 
m > 2 we confirm retrospectively that our solutions possess the correct asymptotic behavior. The 
situation is more delicate in the case of the modes m = 0,1, in which the horizon regularity 



conditions (|77j)— ( 81 ) need to be explicitly imposed in constructing the solutions, as we shall discuss 
in Sec. IVl 

The physical boundary conditions at the poles were given in Eqs. (82) and (83). In our im- 
plementation we make use of so-called "ghost points", which lie in zones just outside the physical 
domain (at 9 < and 9 > vr/2), whose values are artificially set to enforce the correct boundary 
conditions on the physical domain. This allows us to use the same finite-differencing molecule 
everywhere within (the physical part of) the grid. In addition, we halve the computational burden 
by evolving for grid points in the domain < 9 < it/ 2 only, and applying the symmetry condition 
(83) on the equator. 



4- Initial conditions 

In Papers I and II we used very primitive initial data, and relied upon the radiative character of 
the equations to dissipate the "junk radiation" that is initially present. In this spirit, we will start 
with the trivial initial data set = = dtv^ . We will show that this is sufficient for m > 2, 
but that more careful consideration is needed for the non-radiative modes contained in m = 0,1. 



5. Decomposition in tensor spherical harmonics 

The linearized equations governing Lorenz-gauge MPs on Schwarzschild spacetime may be sepa- 
rated using the set of ten tensor spherical harmonics. With this approach, the angular dependence 
of the equations is completely removed, leaving (for each l,m where I > m) a set of ten coupled 
partial differential equations in independent variables t and r. This set may be further decoupled 
into seven equations describing the even-parity part of the perturbation, and three equations de- 
scribing the odd-parity part [lU H6]. Thus far, we have spurned this decomposition, because our 
key motivation is to develop a method that can be generalized to axisymmetric spacetimes like 
Kerr's, where such a separation is not possible. However, we now give this decomposition for two 
reasons: firstly, to make possible comparisons with known results in the literature; and secondly, 
because this decomposition will later facilitate better understanding of the role of non-radiative 
modes in the Schwarzschild case. 

Following Barack and Lousto [17] [Eq. (8) therein], one makes the decomposition 
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Lm i=l 



1U 

W = rEE r)Y^(9, cP; r), (88) 



where the tensorial harmonic basis functions Y^) lm (6 , (p; r) and coefficients a^ 1 are defined in 



Sec. IIB of [T7]. This should be compared with the m-mode expression, Eq. (47). It is straightfor 
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ward to construct the Ira modes 7i|^(t, r) from the m-modes u^g (t,r,0) by integrating over 9. In 
Appendix D we give explicit expression for hi m (t, r) in terms of the m-mode variables. 

IV. NUMERICAL RESULTS: PART I 

In this section we present a selection of numerical results. We show that the numerical simulation 
evolves towards equilibrium for modes m > 2, but that the lower modes (m = and m = 1) are 
plagued by instabilities which grow linearly with t. In Sees. IV] and |VI[ we investigate how these 



instabilities arise, and propose a way forward. In Sec. VII we present some results for the total 
GSF. 



A. m-modes 

Let us consider the results of a typical "run", i.e., a single simulation with a given ro and m, 
and with a particular set of numerical parameters {num.} = { At, Ar*, AO, T rt , Tq, t max }, where 
t = imax is the physical evolution time. We choose to present the results of the run by plotting 
data along three slices (as in Papers I & II): (i) t = t maX i = tt/2, i.e. the equatorial plane, (ii) 
t = tmax ; f = ^0) i-e., from pole to pole and intersecting the worldline, and (iii) r = tq, 9 = ir/2, 
i.e., as a function of time along the worldline. 

Figure [2] shows numerical data for the components of the MP as a function of r* at 9 = it/ 2, 
t = tmax [slice (i)]. Outside the worldtube we plot the real variables 



-( m ) 1 



%'(t,r,0) = Re [u^\t,r,9)e mnt \ , (89) 

and inside we show the punctured version u^ m \ The worldtube is visible as a "trough" around 

ro = 7M (r*o ~ 8.832M). In the asymptotic limits r* — > ±00 the original (complex) m-modes 
go as ~ exp(±imf2r*). This behaviour is apparent as fixed- wavelength oscillations in the numerical 
data for . Figure Jij shows numerical data as a function of polar angle, crossing the worldtube 
at a fixed time t = t max [slice (ii)]. Figure [4] shows numerical data along the worldline as a function 
of time [slice (iii)]. The plot makes it clear that the components of the (regularized) MP on the 
worldline approach stationary values. In other words, the MP is "co-rotating" with the particle, 
and the modes have the expected time dependence. 

B. Validation and tests 

1. Gauge constraint violation 

The numerical solution only represents a physical solution of the linearized Einstein equations 
if the gauge constraints are satisfied. In practice, we do not expect the constraints to be identically 
zero, but rather to be "small" in some sense, and to approach zero as the resol ution i s imp roved. 
Fig. Id shows the profile of gauge constraint violation [defined in Eqs. (A11)-(A14)] as a 



function of r* on the equatorial plane 9 = it/2, in a run with m = 2, ro = 7M at t = 250M. It 
makes sense to compare the magnitude of (which is itself dimensionful) to some norm involving 
the MP's gradient, and in Fig. 5 we have chosen as norm the quantity |<9t/i( m )| = |mO/i( m ^|, where 
/i^" 1 ) is the trace of the m-mode MP. The plots show that the (normalized) constraint violation 
is maximal near the particle's position, but even there it remains quite small: ~ 10 -4 for our 
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FIG. 2: Slice (i): The m — 2 mode of the (trace-reversed) MP as a function of tortoise coordinate r*, for 
9 = 7r/2, r = 7M and t — i max = 300M. The worldtube covers the region r * — r r>> /2 < r% < r * + r r „/2 
where r * ~ 8.8326M and T rt = 8M. The curves show the components u^™ -2 '' [defined in Eq. (89) and 
proportional to the trace-reverse perturbation, Eq. (47)] outside the worldtube, and the punctured version 
" x ' jS> inside the worldtube. The t9, r6 and 9(j) components are not shown as they are zero on the 



equatorial plane by symmetry, Eq. (83). 



Angular profile of metric perturbations 




FIG. 3: Slice (ii): Showing the m = 2 mode of the MP as a function of polar angle 9, for r = vq = 7M, 
^max = 250Af. The worldtube covers the region ir/4 < 9 < Sir/A. The ten curves show the components 
^a"? -2 -* [defined in Eq. (89)] outside the worldtube, and the regularized version m^™ 1-2 ^ inside the worldtube. 



20 




FIG. 4: Slice (iii): Showing the m = 2 mode of the (regularized) MP on the worldlinc r = r = 7M, 9 = ir/2 
as a function of simulation time t. 

highest resolution. More importantly, the constraint violation diminishes as the grid resolution is 
increased. 



2. Decay of junk radiation 

Our trivial (unphysical) initial data create a burst of "junk radiation", which (we find, for 
m > 2) dissipates away with time as the numerical solution approaches a stationary state. An 
example is shown in Fig. [g| The figure shows the real quantity fij- rn ~ 2 ^ defined in Eq. (27) (recall this 



represents the combined contribution from the two modes m = ±2) as a function of t at late time. 
This quantity is constant (non-oscillatory) in the physical solution. We see that junk radiation in 
the numerical solution superposes decaying oscillations of frequency mQ on the physical m-mode 
solution. The amplitude of these oscillations is quite small in the example shown: note the vertical 
scale. [The origin of oscillations in the junk radiation, which is not sourced by the particle, is 
simply the fact that is evaluated in a system co-rotating with the particle — note the factor 



If we were solving the actual perturbed Einstein equations, then we would expect the decay to 
equilibrium to follow the familiar power-law pattern characteristic of generic vacuum perturbations 
of Schwarzschild spacetime at late time: individual radiative multipoles die off as oc H 2i + 3 ), and 
since the lowest radiative mode contained in a generic m-mode perturbation is I = \m\, the \m\- 
mode vacuum MP would be expected to exhibit a oc t~( 2 l m l+ 3 ) decay tail at late time. Indeed, in 
the scalar-field case of papers I and II we have observed a relaxation pattern consistent with this 
rule. Here, however, crucially, we are not solving the perturbed Einstein equations but rather the 



system (52), which is not equivalent to the Einstein equations unless the constraint is satisfied. The 
junk radiation does not generally satisfy the constraint, and thus it does not represent a physical 
solution of the Einstein equations. The rate of decay of these junk solutions is a priori unclear. 
To understand the approach to equilibrium in our system, we experimented with vacuum runs 
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FIG. 5: Gauge constraint violation. The plots show the amplitude of Z^™' [defined in Eqs. (All |— (A14)], 
normalized by \dth^ mS> | = ImQh^ |, where is the trace of the m-mode MP — we denote this normalized 

quantity by Za = Za / \mTlh^ m ^ | . The 3 panels display the t, r and <p components. Each panel shows 
runs at various resolutions n, where Ar* = M/n = At and AO = ir/(8n). The data shown is for the specific 
case m — 2, ro — 7M, t — 250M, and given as a function of r* in the equatorial plane 8 — ir/2. The spurious 
values very near r + — r * should be disregarded; they arise because the retarded field (and hence Z a ) is not 
defined on the worldline. 



starting with some Gaussian (non-Lorenz gauge, non-static) initial data. We found that, at late 
time, the system undergoes power-law decay in which (for m > 2) the BL metric components die 
off as i - ^!" 1 ! -1 ), i.e., significantly slower than the decay of a scalar field in vacuum. We also found 
that certain cancellations occur at late time, so that the gauge constraints Z a and the trace h 
decay as t~( 2 l m l +1 ) — faster than the metric components but still slower than a scalar field. We 
have confirmed the behaviour seen in the 2+ ID simulations by also evolving our gauge-damped 
Z4 system in 1+1D (as in Ref. |17j). Our vacuum 1+1D evolutions for specific multipoles / > 2 
(and m < I) showed t~( 2 ' -1 ) decay tails for the MP, and i~( 2i+1 ) for Z a and h, consistent with our 
2+1D results. 

The above (empirically deduced) relatively slow decay of junk radiation is an unfortunate feature 
of our constraint-damping formulation. It might be possible to improve this feature by controlling 
the form of the constraint-damping terms (oc k) at large r, which is where late-time decay tails form 
through backscattering. We leave this for future study. In the case of circular orbits considered here 
we may reduce the adverse effect of slow junk decay simply by averaging the late-time solutions 
over a complete wave-period T m = 27r/(mf2). Assuming the junk radiation has the late-time form 
~ e im with k as indicated above, such an averaging procedure effectively suppresses the 
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FIG. 6: Power-law decay of junk radiation. The plot shows numerical data for the self-force mode F t f m , for 
to = 2 and vq = 7M, extracted from a time-domain simulation up to i max = 500M [the real quantity F^ m \ 
defined in Eq. (27), is formed by adding together the contributions from the two modes ±|m|). The small 



damped oscillations (note the vertical scale varies by less than 0.04%), with frequency 2fl and an envelope 
scaling as t~ 3 , are due to the power-law decay of junk radiation generated by unphysical initial data. The 
equilibrium value (dashed line) may be estimated by averaging over a complete wave cycle, as described in 
the text. 



amplitude of oscillations (e.g., in the real quantity ) by a factor ~ 2ir 2 [k{k + l)]~ 1 (t/T m ) 2 ~ 
27r 2 (i/T or b) 2 at leading order in t/T m , where T or b = 2ir/Q is the orbital period and where in the 
second expression we have approximated k ~ 2\m\. This is a significant gain for t 3> T or b. We 
have implemented this ad-hoc procedure in our analysis; the outcome is illustrated in Fig. [6| 



3. Convergence of results with grid resolution 

In vacuum, our finite-difference scheme is 4th-order accurate, in the sense that, for grid spacing 
AAt, AAr* and AA#, the discretization error scales approximately as A 4 . In our case, however, the 
irregularity of Sffi e ^ at the worldline disrupts the global convergence of the finite difference scheme. 
As in Paper I, we find that a 2nd-order puncture formulation leads to a global discretization error 
that scales with A 2 In A at leading order. 

In order to improve our estimates of the physical results, we fitted an appropriate model to 
data extracted from simulations at various resolutions. The fitted model allows us to extrapolate 
to zero grid spacing, A — > ("Richardson's deferred approach to the limit", [96]). We applied the 
model 

X(A) = X{0) + q n A 2 In A + c 2 A 2 + c 3 A 3 , (90) 

where c\ n , C2 and C3 are numerical coefficients to be determined, and X stands for u a ^ ^ or ' . An 
example is shown in Fig.[7j for the mode m = 3. We compare the extrapolated value with a highly 
accurate value which was obtained using a frequency-domain calculation by Akcay [35] (summing 
over all frequency modes and all Z > 3 for m = 3). In this typical example, after extrapolation 



with the model above, Eq. (90), we find agreement up to a fractional error of ~ 2.5 x 10' 
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FIG. 7: Convergence of time -do main m-mode data with grid resolution. The left panel shows numeri- 
cal data for p^ m ~^ [cf. Eq. ( 27 )] , extracted on the worldline at ro = 7M and plotted as a function of 
simulation time i, for a range of grid resolutions n, with Ar* = MJn — (5/4) At and A8 = 7r/(8n), 
for n 6 {4,6,8,12,16,20}. The numerical data converges with finite-differencing error oc n~ 2 \nn. The 
right panel illustrates our Richardson extrapolation procedure. Red circles show the values extracted at 
imax = 250M for grid resolutions corresponding to n E {4,6,8, 12, 16,20} from right to left. The blue line 
shows the line of best fit, according to the model given in Eq. (90 1. The black square in the right panel, and 
the straight (blue) dashed line in the left panel, mark the highly accurate value of obtained using a 

frequency-domain calculation |35j . The relative difference between our extrapolated value and the accurate 
frequency-domain value is ~ 2.5 x 10~ 5 . 



4- Comparison with Im-mode results 

Figure M shows numerical data for the angular profile of the (unregularized) retarded MP. More 
precisely, it shows the real and imaginary parts of u ™ e im , for m = 2 at r = vq = 7M. The plots 
show that the components tt, t<j) and <jxj) diverge as the worldline is approached. It can be checked 
that this divergence is logarithmic, i.e. ~ In 1 6* — vr/2|. The other components (for example, tr in 
Fig. [8]) are regular across the worldline. 

We may extract numerical estimates for Im modes from m-mode data by projecting the latter 
onto the tensor spherical harmonics, and using the relations given in Appendix [D] With this 
approach, we have validated our results against previous frequency-domain studies [35]. Table 
[I] shows typical numerical results from this process. Here we have extracted the value of the 
I = m = 2 even-parity (i = 1 ... 7) modes on the worldline at r = ro = 7M. The small discrepancy 
between the m-mode estimate and the (highly-accurate) Im mode results is at the expected level 
of discretization error of our time-domain scheme, ~ 5 x 10 -5 fractionally. 



5. Convergence of the m-mode sum 

As established in Paper I (and in Ref. [20J), with a 2nd-order puncture implementation we 
expect the modes of the residual field u^!) m \ and the modes of the self-force Fr m \ to fall away 
as ~ m -2 in the large-m limit. The dissipative components of the GSF (Ft and Fa, for circular 
orbits), which we expect to be exponentially convergent with m, are exceptions to this rule. The 
individual modes of F t and F^ have a gauge- and puncture-independent interpretation, in terms 
of the rate of loss of energy and angular momentum through gravitational wave emission. On the 
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FIG. 8: The angular profile of the (complex) retarded field variables u^g e' for components a/3 given by 
(clockwise from top left) tt, t<f>, <jxj) and tr. The real parts of the former three components are logarithmically 
divergent at the worldline. Here m = 2 and rp = 7M. 



(0 


1 = 2, m = 2 


1 


3.12437 -0.26317i 
3.12457 -0.26316i 


2 


-0.23121 + 0.97576i 
-0.23121 + 0.97577i 


3 


3.79704 +0.44010i 
3.79727 + 0.44011i 


4 


-0.92491 +9.42906i 
-0.92491 +9.42918i 


■5 


-2.33102 -2.52793i 
-2.33103 -2.52790i 


G 


1.54687 +0.60066i 
1.54684 + 0.60065i 


7 


-5.33209 -5.21910* 
-5.33189 -5.21903i 



TABLE I: Example of numerical validation of m-mode results against Zm-mode results for the case I = m = 2 

~ (*) J I 

7M. The index i labels individual tensor-harmonic modes h lm (r = ro) [cf. Eq. (88)], with 



and ro = 

i = 1, . . . ,7 corresponding to modes of even parity. In each entry, the upper value is extracted from our m- 
mode scheme, using Eqs. (Dl |-(D10). The lower value is obtained from a highly accurate frequency-domain 
calculation |35j . The small discrepancy is consistent with the expected level of error in the m-mode data. 
Similar accuracy is found for odd-parity modes (e.g. / = 3, m = 2, i = 8 . . . 10) and for higher modes, m > 2. 
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FIG. 9: Convergence of the m-modc sum. The left panel illustrates, on a log-log scale, the m 2 power-law 
fall-off (at large m) of the m-mode contributions F^ for a circular orbit at r$ = 6M. The dashed line is a 
reference ~ 0.0052/m 2 . The right panel shows, on a semilog scale, the exponential decay exhibited by the 
m-modes of the dissipative GSF component. The dashed line is a reference 0.003 x e~ 1Am . 



other hand, the individual modes of F r are somewhat arbitrary, since they depend on the analytic 
extension of the puncture function away from the worldline (which is implementation-dependent). 

Figure [9] confirms, in the example of ro = 7M, that our expectations for the large- m. behaviour 
of F^ and F^ are met. 



C. The modes m = and m = 1 

In the preceding sections we have shown that, for the modes m > 2, the MP settles into an 
equilibrium configuration, and initial junk radiation dissipates with time. Unfortunately, this is 
not the case for the modes m = and m = 1. 

In earlier 1+1D studies of the GSF, the modes I = 1 and I = could not be obtained through 
the time-domain approach, and were obtained separately using frequency-domain analyses |97| . 
Indeed, a fully-time-domain Lorenz-gauge scheme has not yet been successfully developed. Barack 
and Sago [98] noted that "experimentation suggests to us that the monopole and dipole cannot 
be evolved stably using this [time-domain] scheme. A naive application of the evolution scheme 
yields exponentially growing solutions, and, since our scheme gives us no handle on the boundary 
conditions, the occurrence of these unphysical solutions is difficult to control." 



With our choice of gauge-constraint damping (Sec. IIIB 1) and a Cauchy (t, r*) scheme we do 
not observe exponentially growing solutions. On the other hand, we do observe solutions which 
grow linearly with time t. Numerical data from typical evolutions of the m = and m = 1 
modes are shown in Fig. [TUJ In both cases, it seems that the instability can be attributed to a 
homogeneous (vacuum) mode, which is generically excited by our unphysical initial data. Such 
homogeneous modes appear in the co-rotating variables u^g as oscillations with frequency mO. 



Figure 11 shows the radial profile of a typical evolution of the m = mode, at late times. The 
plots indicate that the components of the MP that grow linearly in t also scale as ~ r* near the 
horizon, and we find that the near-horizon scaling is in fact proportional to the advanced time 
coordinate v = t + r*. In the near horizon limit we find hu ~ ht r ~ h rr ~ —\hog ~ —\h^ = Av 
for some amplitude A which depends on initial data. In other words, these troublesome modes 



demonstrate "ingoing" behaviour at % + [i.e., (dt — d rt ,)h a p — > 0]. Figure 11 c) shows the com 
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FIG. 10: The non-zero components of the m = and to = 1 modes of the (regularized) MP variables 
on the particle's worldline [cf. Eq. (89)] at r = 7M, as a function of time t. Left panel: the m = mode, 
showing linear growth in the tt, tr, rr, 99 and </>(/> components. The t<j) component approaches a stationary 
value. Right panel: the m = 1 mode, evaluated on the worldline (at <fi — Qt). The linear-in-t growth appears 
on the worldline as a growing oscillation of frequency il in components tt, rr, r6, r<p, 99 and cjxp. For both 
to = and to = 1 the trace approaches a stationary value (inset), which suggests that the linear-in-t part is 
traceless. 



ponents of the Lorenz-gauge violation Z a . We observe that the gauge violation remains small, 
suggesting that the growing modes are indeed valid Lorenz-gauge solutions. 

By examining snapshots of the pole-to-pole angular profile, we infer that only the low multipoles 
I < 2 are implicated in the instability. This suspicion is further confirmed by evolving the equivalent 
Im field equations in 1 + ID [Eqs. (91 )-p4|) below]. 

In the next section, we derive explicit analytic expressions for the above linear-in-t modes that 
are disrupting our naive implementation of the time-domain scheme. We show that these modes 
are pure (Lorenz)-gauge modes, which are closely related to "scalar-type" (traceless) gauge modes. 
In Sec. |VI| we then address the challenge of "stabilizing" the m = and m = 1 evolutions, and in 



Sec. VII we display some final results for the total GSF and compare with the literature. 



V. LOW MULTIPOLES: ANALYTIC CONSIDERATIONS 



In this section we consider the non-radiative multipoles of the MP, i.e., I = and I = 1 in 

Physical solutions for these modes have been discussed before, 
[17] [521 [57] [99] . Here, we give a further analysis which focuses on two issues: 



the decomposition of Sec. |III E 5 
for example in Ref. 

(i) the presence of conserved mass-energy and angular momentum in the low multipoles of the 
MP, associated with Killing vectors of the background spacetime (see the covariant derivation in 
Sec. |II E| ), and (ii) the existence of linear-in-t (and Lorenz-gauge) gauge modes, which are regular 
on 7i + (in the sense of Sec. Ill D ) as well as at infinity. A key point is that such gauge modes, which, 
at any finite time, are well-behaved everywhere, cannot be eliminated using boundary conditions 
alone. To eliminate these modes we must impose additional conditions; for example, the condition 
that the monopole perturbation is static (i.e. dth^ u = and h tr = 0). 



Some parts of the analysis below, especially from Sec. VA3 onwards, will take a more general 
form (less tied to the 1+1 D context), with the idea of preparing the ground for a similar analysis 
in Kerr. 
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FIG. 11: Features of the linear-in-i instability in the m = mode, illustrated here in the equatorial plane 
at t = i max = 250M for circular orbit at r$ = 6M. (a) Components of the MP u^~°\ showing that in the 
near-horizon regime (r* < — 10M) the components tt, tr, rr, 89 and cfxfi scale linearly with r*. (b) The time 
derivative of the components, showing that in the near-horizon regime the same components scale linearly 
with t as well; in the text we conclude that the growing mode is "purely ingoing" at the event horizon, 
(c) The (normalized) gauge-constraint violation, which appears to be small and w h ich d iminishes as grid 
resolution is improved. Here we show the amplitude of Z^™' [defined in Eq. ( All )-( A14)], normalized by 



|jlfft(m=o)| w h ere h(m=o) is the trace of the 



MP. 



A. Monopole perturbation 

The I = m = monopole mode is spherically symmetric and of even parity; physically, it 
describes the mass perturbation due to the particle (up to a gauge). In the Lorenz gauge, this 
mode is governed by four field equations and two gauge constraints. Analytic solutions for the 
monopole were previously obtained in Sec. Ill of and in Sec. HID of [T7] for circular orbits. 
Here we present an alternative, but closely related, analysis. 



1. Monopole equations 

Let us begin by stating the monopole equations in terms of the spherically-symmetric part of the 
relevant m-mode variables, i.e., u a /s = \ Jq sin Gu^~°\t, r, 9)d9. The nonvanishing components 

— (i) 

are related to the Im modes h lm of Eq. 

and uqq = = ah^, where a = (16-7r) -1 / 2 . The four Lorenz-gauge equations without constraint 



88) 



via u tt = a{h\^ + h\^ 



), U tr 



1,(2) 
an lm ' U rr 



h {3) ) 



28 



damping are 



D utt + 



2{u u - u rr ) , 4(u 4r - u' tt ) 4/(u« - U rr ) 4:f 2 U d9 _ 4/i£ 



+ 



+ 



+ 



To 



fo$(r* - r *), (91) 
(92) 



2 2{u tt + Urr ~ 2u' tr ) 2f 2 U tr 
U Utr + 7. 7, — = U, 

2 2(u rr - u tt ) 4(tt tr - u' rr ) £f{u tt - u rr ) 4/ 2 u rr 4/ 2 n e g 4/ 3 ngg _ 



D 2 u ee + 



2(u tt ~ U rr ) 2f(u rr + Ugg) _ 4/ 2 ttgg 



2u£ 

(ro ~ 2) 



r" 



2 fo$( r * ~ro*), 



r- 



(94) 



where D 2 = —d 2 + d 2 t — 2//r 3 , and u and u' denote differentiation with respect to t and r*, 
respectively. Recall ro is the orbital radius and /o = 1 — 2M/rQ. The monopole gauge constraints 
may be written as 



_ Kr - iht) fu tr n 
<^t = 1 — — U, 



r 



{u' rr ~ u tr) (utt ~ U rr ) fu rr _ 2f 2 U d _ 



(95) 
(96) 



To obtain the dynamical equations with gauge constraint damping of the form given in Sec. IIIB 1 



one then takes the combinations Eq. (|9l|) + 4Z t /r, Eq. (j92J) + 2(Z t + Z r )/r, Eq. ^3j + AZ r /r and 
Eq. ([941). 



U. Mass-energy condition 



In Sec. 



HE 



we obtained a relation between the particle's energy u£[= —Q(X?L)] and a closed 
2-surface integral T over an antisymmetric tensor F a p constructed from the MP, its derivatives, 
and the background timelike Killing vector X%s [Eq. (28)]. We now specialize to a Schwrazschild 



background and to a circular orbit of radius r = ro, choose our 2-surface to be a 2-sphere of r=const 
(and i=const), and define the "energy" functional E[u a ^;r] = —T{X? t ydll r ). We then have 



E\u, 



Q/3\ 



F tr dn 



u£, r > r , 
0, r < ro. 



(97) 



An explicit expression for E[u a g] is obtained by substituting for F tr from Eq. (28) with X a = 5", 
and noting that the surface integral picks out the monopole part of the MP: 



E[u a8 ] 



■\rr x 



u 



f 



tt - u tr ~ -U tt 



2(u tt 



It* 



f(u tt + U rr ) 3(u tt ~ U rr ) , 2f 2 U 9e 



+ 



(98) 



where, in going from the first line to the second, we have used the Lorenz-gauge constraint (96) 



to eliminate iit r . It is straightforward to use the vacuum version of Eqs. (91)-(96) to confirm 



explicitly that dtE[u a g] = = d r E[u a p], i.e., that the MP combination on the right-hand side of 
(98) is indeed constant for r ^ ro, as expected. 



Eq. (97) gives a necessary condition for the MP u a g to represent a physical monopole solution. 



The monopole piece of our m = numerical solution must satisfy this condition. 
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3. Static homogeneous Lorenz-gauge solutions 



Now let us consider the family of solutions to the homogeneous part of the monopole equations 
(|9l]H96). With the staticity condition ut r = 0, we have 3 independent second-order field equations 
[Eqs. (91), (93) and (94)], and a single nontrivial gauge condition [Eq. (96)]. Since the gauge 
condition (96) gives uqq algebraically in terms of u tt and u rr (and the first derivative u' rr ), the 
system effectively reduces to a set of two second-order ordinary differential equations. Hence we 
expect the complete basis of homogeneous solutions to be 4-dimensional. Below we will construct 
a complete basis of four homogeneous solutions. We will characterize each basis solution by its 
mass content, by whether it is traceless or tracefull, and by whether or not it is regular at the 
future horizon and at spatial infinity. Of our four solutions, only one will possess a non-zero 
mass-energy; the rest will be pure-gauge solutions. For notational simplicity, in the rest of this 
section we usually adopt the convention M = 1. 

Solution A: — Our first solution, naturally called the "conformal" solution, is given by 



11 [IV 



(99) 



This is clearly a valid Lorenz-gauge solution because the covariant derivative of g^ u is zero (hence 



9iiu' u = and also Og^ = 0), and the background is Ricci-flat (hence 2R\ v g\ a = 2R fll/ = 0). 
This solution has a (constant) nonzero trace (M A ) = 4), and a nonzero mass energy E[u a p\. It is 



regular at the horizon [in the sense of Eqs. (77)-(79)], but note that it does not fall off to zero at 
infinity (but approaches a constant value there). 

Solutions B and C: — Next we consider the "scalar" pure-gauge Lorenz-gauge solutions, given 
by h a p = £a;/3 + £|8;af) where £ a = $ ;a for some scalar field <£. Note the MP trace is h = 2D$. The 
Lorenz-gauge condition, together with Ricci-flatness, implies that 



h, 



2D£ Q = 0, 



(100) 



hence such solutions have a constant trace. The ansatz $ = $(r) leads to two independent scalar 
pure-gauge solutions given by 



Kr 3 -8)/(fr 2 ), 



(101) 



and 

(with all other components vanishing), of which the first is tracefull (h^ 



traceless (h^ = 0). The corresponding MPs are given explicitly in Eqs. (106) and (107) below. 



(102) 
6) an d the second 



These solutions are massless. We will consider their regularity properties in the next subsection. 

Solution D: — Finally, let us consider a more general pure-gauge homogeneous solution, derived 
from a gauge displacement 



(103) 



which is (up to a multiplicative constant) the most general form of £ Q consistent with a static MP. 
The Lorenz-gauge condition gives an inhomogeneous second-order ordinary differential equation for 
£m)( r )- The above scalar modes £r and & are two independent solutions of the homogeneous 
part of this equation, and as a particular solution to the inhomegenous equation we take 



(D) 



3r 2 



[(r 2 + r + 4)r 



ln/ + 81nr] 



(104) 
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The corresponding MP in given in Eq. (108) below. This solution is tracefull {hS D ' = 2 In/) and 
massless. We will consider its regularity properties in the next subsection. 

To summarize, we have constructed a complete basis of independent homogeneous solutions 
A-D to the static monopole equations. Writing h = {h, E,h t t,h rr ,r~ 2 hQQ = (r sinfl) -2 /^} for 
any of the solutions A-D, we have, explicitly, 

h A = ^{4,l/2,-/,r\l}, (105) 

ti B = n{6,0,-2fP(r)/r 3 ,2f- 1 Q(r)/r 3 ,2fP(r)/r 2 } J (106) 

tic = M0,0,-2/r 4 ,-2(2r-3)r 2 A 4 ,2/r 3 }, (107) 

ti D = fi{2ln f, 0,- ^[rW(r)+rfP{r) In f -81nr], 

2 



3/ 2 r 4 

_2 
3r 



[-rK(r) + L(r) In / + 8(2r - 3) In r] , 



2 i r (p( r )- r ) + ( r 3 - 8 )ln/-81nr]}. (108) 



In these expressions 



P(r) = r 2 + 2r + 4, (109) 

Q(r) = r 3 -r 2 - 2r + 12, (110) 

W(r) = 3r 3 - 7r 2 - r - 4, (111) 

K(r) = r 3 - 5r 2 - 5r + 12, (112) 

L(r) = r 4 - 3r 3 + 16r - 24. (113) 



4-. Static inhomogeneous Lorenz-gauge solution 



We now attempt to construct a unique physical monopole solution, fulfilling the following list 
of requirements, (i) It is a solution to the inhomogeneous equations (91)-(94). (ii) It is a Lorenz- 
gauge solutions, i.e., Eqs. (95) and (96) are satisfied, (hi) The MP is static in the sense that 
dth a R = and hu 



for 



(iv) The MP is continuous across 



The mass- 



energy condition (97) is satisfied, (vi) The MP is regular across T~L + (Sec. HID), (vii) The MP is 



regular at r — > oo, in the sense that h^y/^^y (no summation) falls off to zero at least as ~ 1/r. We 
will find (in confirmation of [T71 [211 EZ] ) that it is not possible to construct a solution that meets 
all such requirements. 



Let us first consider horizon regularity. It is easy to see, using the the conditions of Eqs. (77) 



(79), that solutions A and B are regular at the horizon, whereas C and D are not. However, the 



conformal solution A has mass-energy (and it is the only solution with this property), so by above 
requirement (v) is it not allowed within the interior of the orbit. Thus we may immediately write 



Lor 



/lint = a B h B , r < r , 

ti ext = 2£hA + b B h B + b c h c + b D h D , r > r , 



(114) 



with four coefficients a B , b B , be and b n to be determined. Here the coefficient of hA in the exterior 
is determined by the mass condition (97), recalling, again, that Ha is the only massfull solution, 
with E = fi/2. The continuity requirement (iv) imposes 3 additional conditions, giving 



a B - b B 

be 
b D 



-|[(r -3)ln/o-4ro/o] 

-| [r 2 + 20r - 64 + 8(r - 3) mr ] 
-a(r - 3), 



(115) 

(116) 
(117) 
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where fo = 1 — 2/ro, and 



a=-*—. (118) 
rofa 

Now let us consider the asymptotic behaviour at r — > oo. We find 

2/xa + 0(r- 1 ), (119) 

fcrr, r" 2 ^, (rsinfl)" 2 /^ ~ 2 ^ a ( 4 ^ ~ 9 ) + ^ + Q{r -iy (120) 

The leading (constant) term in the spatial components (rr ,98 ,(f>(f)) can be removed by fixing the 
remaining degree of freedom to be 

b B = -a(4r - 9)/3, a B = a[l - (r - 3) ln/ ]/3. (121) 

However, the it component of the perturbation remains asymptotically nonzero, which cannot be 
remedied without violating at least one of the other requirements (i)-(vi). 

It is straightforward to verify that the static solution given above is identical to the solution 
presented in Sec. HID of Ref. [T7] (or the solution inferred from |52j), as expected. Here, for the 
first time, we derived this solution from a complete basis of (static) solutions to the Lorenz-gauge 
monopole equations. Having such a complete basis at hand is useful in general, because it allows 
the construction of a physical mass-perturbation solution also for non-circular orbits. [Indeed, the 
above complete basis was previously used by one of the authors, in [98], to construct a physical 
monopole solution for generic bound (eccentric) orbits; however, Ref. |98j did not give the complete 
basis explicitly, as we do here.] Note also that our construction avoids the explicit imposition of 
jump conditions for the MP derivatives at ro, but it can be checked that these are satisfied. In 
our construction we have replaced the jump conditions with a "mass condition", with the foresight 
that in the Kerr problem one can still impose the latter, but one cannot impose explicit monopole 
junction conditions (because the "monopole" piece cannot be separated out in Kerr). Indeed, many 
of the details of the above analysis are transferable to the Kerr case, as we shall show in paper IV. 

To obtain an asymptotically regular mass-perturbation solution, it is most natural to relax the 
Lorenz-gauge condition. This approach was taken in Ref. [57], via the introduction of the simple 
non-Lorenz gauge transformation 



with = —fj,atSf, so that 



*C L = C+C> ( 122 ) 

5hf t L = 2fiaf, (123) 



and the other components are zero. Recalling Eq. (119), we see that the new perturbation + 
Shf t L is C(r _1 ) at r — t- oo as desired. We note, however, that this new perturbation now fails to be 



regular at % + : From Eq. (79) we see that its Eddington-Finkelstein RR component diverges there 



as ~ . An alternative choice, £^l = ~f ia ('t + r* — r)S^, leads to 

— NL — NL 

5h tt = 2fj,af, 8h tr = 2/iaM/r, (124) 

with all other components unchanged. This MP is regular at T~L + , but it is not static in the sense 
of condition (v). We have not been able to find a simple transformation away from Lorenz gauge 
leading to a monopole solution that is static, continuous, globally regular, and has the correct 
mass-energy content. 
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Staying within the Lorenz gauge, it is in fact possible to write down a Lorenz-gauge monopole 
solution that is static, continuous and globally regular - but has the wrong mass-energy E. This 
solution, first obtained by Berndtson 100 . is given by 

fy* no = h^ no -2ah A + ah B . (125) 

While this solution is nonphysical [because it corresponds to a mass perturbation from a non- 
geodesic particle with energy //£ — 2/ia = fx£ (rp— 3M)/(ro— 2Af)], it isauseful solution nonetheless, 
because it can arise in time-domain simulations (which usually seek globally regular solutions). If 
we identify this nonphysical solution in the numerical data, we can easily "correct" the mass-energy 



by applying Eq. (125) in reverse. We will come back to this point later. For future reference, we 



shall refer to /i^° r no as the "asymptotically-flat" Lorenz-gauge solution. 



5. Nonstatic and nonstationary (linearly growing) monopole solutions 



In Sec. IV C we saw that time-domain evolutions of the Lorenz-gauge monopole equations can 
excite pure-gauge modes that grow linearly with time t. In what follows we give an analytic 
description of such modes. For notational brevity we set fj, = 1 = M throughout this subsection. 

First, observe that if the above staticity condition (v) is relaxed, then the monopole MP is no 
longer unique, even if all other conditions are satisfied. Consider, for example, the Lorenz-gauge 
scalar perturbation $ nonstat = |i + ln/, for which £™ nstat = ^nonstat = 2/(r 2 /), 0, 0] . The 
nonzero components of the corresponding MP are 



nonstat 



nonstat 



nonstat 



h 



r 2f> 

4(2r - 3) 
r 4 / 2 ' 



^r tat = h H /sm 2 e=^. (126) 



This homogeneous, Lorenz-gauge MP is regular on Tl + according to criteria (77)— (79), and also 

.non 



asymptotically flat. It is stationary in the sense that <9(/i^° nstat = 0, but it is nonstatic since 



ht r 7^ 0. This MP is traceless, /j nonstat = 0. As we demonstrate below in Fig. 12, this nonstatic 
mode features in our 2+1D numerical m = solutions. 

We seek a solution that satisfies all criteria (i)-(vii), except staticity, and for which the metric 
components grow linearly in time. We are guided by the observation that our linear-in-t numerical 
solutions seemed to have perfectly stationary (t-independent) mass-energy E and trace h. This 
strongly suggests that the mode implicated in the linear growth is pure-gauge and of a constant 
trace (or traceless). We will now construct analytically such a linearly-growing Lorenz gauge mode, 
which is globally regular (at any finite time). 

Let us first consider the scalar gauge mode $ = |fln/, which is a simple monopole solution of 
□<& = 0, and hence it generates a traceless Lorenz-gauge solution. This leads to h aa = t x h^ a (no 



sum), where h^ a is the homogeneous solution C given in Eq. (107), and ht r = (2 — ln/)/(r 2 /). The 
diagonal components of this solution grow linearly in time, whereas the tr component is stationary. 
Now consider a second scalar mode, constructed from <I> = t 2 + R(r) where R(r) is chosen so that 



= const. This generates a MP of a constant trace [recall Eq. (100)], which has stationary 
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diagonal components, and a tr component that grows linearly in t. By taking a certain linear 
combination of these two modes we construct the gauge vector £jj n = (£| m , £j! n > 0, 0), with 



^ in = ln(2/)+t/2 + 



13 



flin 



2t 

r 2f 



+ 



+ 3r 2 + 12r + 241n(/r) 
6r 2 / 



6/' 



(127) 



This generates a homogeneous Lorenz-gauge monopole MP hr^i, whose nonzero components read 



Jin 



-r 4 + At + r 2 + 4r + 8 m(r/) 



• lin _ *+i + 21n(2/) 
ft* 



''tr 



h 



lin 



r 2 / 

4i(2r - 3) + 5r 2 - 12r + 8(2r - 3) ln(r/) 

r 4 ' 



-2, lin 4t + r 2 + 4r + 81n(r/) . 2 

r = 5 = (rsm0) z h 4 



(128) 



We note that the origin of time is essentially arbitrary here, so t — > t + c also yields a solution, 



which is equivalent to adding a multiple of the stationary traceless solution (126). 



It can be checked that the pure-gauge solution (128) is regular at according to the criteria of 
Sec. HID It is also regular at r — > oo in all components except tt, for which /i^ n ~ 1 + C(r _1 ). It is 
a pure-gauge solution with zero mass-energy, E = 0, and a constant trace h = 1. By combining 



with the inhomogeneous static Lorenz-gauge solution hh?* no of Sec. V A 4, it is possible to construct 



a Lorenz-gauge monopole solution with the correct mass-energy, which is both horizon-regular and 
asymptotically flat. The metric components of this solution grow linearly with t, but its trace is 
stationary. These features are all in common with those of the linearly-growing mode observed 
in our m = numerical simulations. Hence, the solution h^a a good candidate for being the one 
implicated in the linear instability. 

The above suspicion is indeed confirmed by examining the numerical data. An example is shown 



in Fig. 12 In Sec. VI we shall take advantage of the analytic insight thus obtained into the origin 
of the linear instability, in order to find a cure to it. 



B. Dipole perturbation 

A particle in a circular equatorial orbit generates a dipolar / = 1 MP with two parts: an odd- 
parity perturbation in the axisymmetric I = l,m = mode, and an even-parity perturbation in 
the I = 1, m = ±1 modes. The former part is attributed to the angular momentum perturbation 
due to the particle, and the latter is associated with the recoil motion of the large black hole at 
O(n). We consider each of these two parts in turn. 



1. Odd-parity mode (I = l,m = 0) 

The odd-parity axisymmetric Lorenz-gauge dipole perturbation is governed by a single dynam- 
ical equation, with two independent static homogeneous vacuum solutions given (up to multi- 
plicative constants) by [T7] h t( p = fir 2 sin 2 9 and h t( j, = — 2/ir _1 sin 2 # (other components vanish). 
The former solution is pure-gauge, and is generated by the vector £ a = [0, 0, 0, t] (which satis- 
fies the Lorenz gauge condition, = 0). The latter solution contains "angular momentum" 
L[u a /3;r] = — .F(X?U , <9£ r ) = 1. Note that the latter solution may be obtained through a direct 



34 




1 




/ 




/ 




/ 








1 / 

1 / 




1/ 

V 




I 








1 







50 



100 



150 
t/M 



200 



250 



300 



\ 



i A 



V 



\ / 

i / 

150 
t/M 



FIG. 12: Removal of the linear-in-t gauge mode from the m = numerical data. The implicated mode 
is a certain homogeneous Lorenz-gauge mode of a monopolar angular dependence, which we identified 
analytically in Eq. (1281. Here we remove the problematic mode (for the case shown in Fig. [Tl] with 



7'o = 6M) by subtracting a suitable amplitude of this analytic solution. We fix the amplitude (and "origin 
of time") of the gauge mode to be subtracted by demanding that ht r = ht r = on the worldline (after 
subtraction). Above, the upper plots show the monopole mode for r = 6M, obtained from m = data by 



projection [via Eqs. (Dl )-(D10)], after the removal is attempted. The left plots show the radial profile, and 
the right plots show the solution on the worldline as function of time. The data shows that the "cleaned" 
monopole solution is both stationary (dtu a p) and static (ut r = 0), up to numerical error. The lower 

plots show the difference between the numerical data and the analytic monopole solution h^ no [Eq. (125)] 



obtained by Berndtson |100j which, as argued in Sec. VA4 is asymptotically-regular but has incorrect 
mass-energy. The difference between the numerical and analytical solutions is small, and diminishes as grid 
resolution is improved. Nevertheless, the numerical accuracy remains unsatisfactory and so in Sec. VI we 



develop alternative methods to obtain stable evolutions (Fig. 14 ) 



linear variation of the Kerr metric with respect to the background angular momentum aM (with 
fixed Boyer-Lindquist coordinates), at the limit a — > 0. (This procedure yields a Lorenz-gauge 
perturbation only in the a — > limit.) 



Following the argument of Sec. HE, the physical solution should have an angular momentum 



of L = /j,C for r > ro, and of L = for r < r$. Therefore, 



(l=l,m=U) _ I ^f"~v' /'0J O111 u i ' ^ 'U> (129) 
t( t> ~ 1 _o„/v.-i ■ 9 " 



i (i=l,m=o) _ J -2/i£(r 2 /?o) sin 2 9, r < r , 
-2/i£r _1 sin 2 6, r > vq. 

This matches the solution originally obtained by Zerilli [32] (and reproduced in Refs. \17\ \ 

Within the m-mode scheme in 2+1D we do not decompose into /-modes, and we start with 
trivial initial data. Then it is not guaranteed, a priori, that the numerical solution would evolve 
towards one that has the correct angular momentum content. However, by computing the surface 
integral L\u a p]r\ we can "measure" the angular momentum in the numerical m = solution, 
and, if necessary, simply add a multiple of the homogeneous solution hts = — 2r _1 sin 2 9 (which is 
regular on in order to "adjust" the angular momentum. In practice, we indeed find in our 

implementation that the m = numerical solution does not carry the correct angular momentum, 
and we apply the above simple cure. 
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The odd-parity dipole does not contribute to the linear instability observed in our m = 
solutions. This is clear from the fact that the unstable mode has a monopolar (^-independent) 
angular profile. 



2. Even-parity modes (I — \,m = ±.1) 



The time-domain simulations of Sec. IV C (see Fig. 10) suggest that the I = 1, m = ±1 



even-parity modes are susceptible to a linear gauge instability. Let us explore the origin of this 



instability. As in Sec. VA5, let us begin by considering the static traceless scalar (Lorenz-)gauge 
mode £ a = with □$ = 0. This admits separable I = 1, m = ±1 static solutions of the form 
= ij)(r)y±, where y± = smBe^ 1 ^ and hereafter the label ± corresponds to m = ±1. Two 
independent radial solutions are 

-0! = r -M, and ip 2 = 2M + (r - M) In /. (130) 

Now let us consider a solution of the form &± = t i ip{y±. This is also a traceless scalar Lorenz- 
gauge solution (satisfying □<£ = 0), and it gives rise to a linear- in-t perturbation. This perturbation 
is not regular on T~L + , and nor is it regular as r — > oo. However, we can construct a globally 
regular, time-growing solution by combining it with the static scalar mode ip2y±, and also with an 
"electric" -type gauge mode ^ oc 5*y± (we use here the language of Ref. [99]). The result is 

£ = -2(r-2)<Sp; ± + *±, where $± = [t(r - 1) + 2 (2 + (r - 1) In /)] y±. (131) 

This gauge vector gives rise to a Lorenz-gauge MP which is globally regular and whose components 
tt, rr, r9, rcf>, 96 and <f><j> grow linearly in time, although its trace remains static, h = 0. It therefore 
has all the characteristics observed in the linearly-growing mode of the numerical m = 1 solution. 



In Fig. 13 we show that the linearly-growing part of the numerical data (obtained by applying a 



filter) is indeed that given by the analytic solution in Eq. (131). 

Note that, since the problematic growing gauge mode satisfies physical boundary conditions at 
any finite time, it cannot be eliminated using boundary conditions alone. Regularity at T~L + also 
means that schemes employing alternative time slicing (e.g., hyperboloidal |101|. 1102] or double- null 
constructions) would not eliminate such instabilities, although these routes remain to be explored. 

VI. LOW MULTIPOLES: GENERALIZED LORENZ GAUGE AND STABILIZATION 



In Sec. II C| we discussed the formulation of the linearized Einstein equations in a generalized 



Lorenz gauge (GLG), which is defined by a gauge driver function (= in the Lorenz-gauge case). 
In this section, we make use of the freedom to choose the gauge driver to seek a formulation 
which is stable, in the sense that time-domain evolutions of generic initial data are free from linearly 
(or exponentially) growing modes, and which lead to stationary values for u£g on the worldline. 
We then describe a method for recovering the Lorenz gauge solution from the numerical solution 
in the generalized gauge. 

A. The in — perturbation 



The gauge-mode instability seen in numerical results for m = 0, 1 (Sec. IV C) arises even in 
vacuum, for generic (non-Lorenz-gauge) initial data. In the previous section, we argued that the 
instability is entirely in the monopole and even-parity dipole sectors. We may therefore use the 
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Time-derivative of linearly-growing mode in even-parity dipole I = m = 1 




FIG. 13: Numerical check on the relevance of the analytic solution for the linear-in-i Lorenz-gauge mode in 
the even-parity dipole sector. In Sec. VB2 we asserted that the numerical solution obtained in the m = 1 



case is a superposition of the physical solution, with time dependence exp(— iClt), and a pure Lorenz-gauge 
dipole (/ = 1) mode, with linear-in-i behaviour [Eq. ( |131 )]. The linear growth in the numerical solution 

may be estimated by applying a simple filter. The plot above shows Re ^[1 — iQ _1 t?{] c^/i]^ ^ for i = 1 ... 6, 
where h 

very simple linear-in-i profile given by 2hW 



/m denotes the radial variables used in Ref. [17J, and r = 6M. The analytic solution (j 131 1) has a 



-M 5 ) = 2//i< 6 ) oc f/r with fcW = hW = /i (4) = 

where here hfi = d t h^. The plot above shows that the time-derivative of the linearly-growing part of the 
numerical solution is in excellent agreement with the analytic solution (with the difference between the two 
shown in the inset). 



simpler framework of the vacuum 1 + 1D monopole equations [formulated in Eqs. (91)-(96), and 
with suitable GLG corrections] to develop and test ideas for stabilizing the evolution. Once this is 
achieved, we will implement the same ideas in our m-mode scheme. 

First, we have tried out various choices of a generalized Lorenz gauge driver H^, with the aim of 
eliminating the linear-in-i modes without exciting other (e.g., exponentially-growing) instabilities. 
The goal was to obtain a stationary monopole MP that is regular at T~L + and approaches = 
(i.e., Lorenz gauge) at late time. After some experimentation, we focused on the following general 
class of gauge drivers, which gave good results: 

H lx = xn^ 1 where X=~U. (132) 

Here = [1, / , 0, 0] is an ingoing null vector, U is chosen to be one of the metric components 
u g K ( r =0) '}, A is a constant, and re is a positive integer. We settled on the specific 

,(m=0) , _ 1 l A iv-i — Q TV,„ TT - „,( m =°) 



choice U = u tr , A = —1/4 and n = 3. The choice U = u tr guarantees that, if the eventual 
monopole solution is static (i.e., °^ = up to numerical error), then the solution will also be 
in Lorenz gauge, = = 0. 

With the above choice of a GLG, we established experimentally that the sourced 1+1D monopole 



equations also evolve stably, and that the solution is regular at the 



lorizon. In the next step we 
specified as initial data the analytic solution /i^° r no given in Eq. (125), which, recall, is a globally- 



regular Lorenz-gauge monopole (which, however, does not have the correct mass). With these 
initial data, we found that, as the numerical resolution increases, the system approaches staticity, 
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i.e., dtu a /3 — > and Ut r — > 0. As discussed in Sec. V A4 it is then easy to "correct the mass" of 
the monopole solution, by applying the the analytic transformation (125) in reverse. 

With the above preliminary tests passed, the next challenge was to implement the GLG approach 
in the 2+1D setting. Moving to a GLG leads to modifications of the constraint-damped Z4 system, 



Eq. (17). Specifically, a new term B a p is introduced [see Eq. (15)], the constraint damping terms 
C a/ 3 are changed, and the effective source may be modified. However, with the simple choice 
of H a oc utr above, the required modifications are rather simple, because the tr component of 
the m = mode of the puncture function happens to vanish (for circular orbits). This implies 
that the effective source is not modified, and we only have to include extra terms on the left- 



hand side of the evolution equations. The extra terms are generated by B a s in Eq. (15), and 



C a /3 — > Ca/3 — K(n a H/3 + npH a ). At the level of the m-mode equations we evolve Eq. (84) with 



AMi?- 0) = -fr U + X ' + ^ (1 - 3/r)) + C, (133) 
AM^ =0) = -fr(x + x'- 2 4)+C, (134) 



AM%= 0) = AM!r 0) = -fdex, (135) 
AM r T =0) = ~fr (x + x'-y(l- 1/r)) + C, (136) 
AM^= 0) = AM^= 0) = -r( X ~ x'), (137) 

and C = —Afx/r. All other components of AA4^l~ °^ are zero. 

Figure [14| shows sample numerical results for the m = mode, using the scheme outlined above. 
The evolution starts from analytic initial data given by the combination of the monopole solution 



(125) and the dipole solution (129). Initially, the tr component is zero everywhere, but, for a finite 



resolution, it evolves to a stationary solution, which is small but non-zero. Figure 15 shows that, 
as the grid resolution is improved, the stationary value of ut r converges to zero. Hence we are able 

to recover the static m = MP in Lorenz gauge, although with a monopole part h^no having the 



"wrong mass". This, however, is easily rectified by applying the the analytic transformation (125) 
in reverse. 

The above demonstrates that, when appropriate initial data is known analytically, the correct 
m = Lorenz-gauge solution can be obtained through a time-domain evolution. But a natural 
question arises: how can we proceed in the Kerr case, where monopole initial data are not available 
in the Lorenz-gauge, and where the equations do not separate into 1+1D? In other words, if we 
are completely ignorant of the correct initial data, can we still devise a scheme that will evolve 
towards the static (ht r = 0) Lorenz-gauge solution? 

We tested a simple but general approach. First, we evolved the sourced equations starting 
with trivial (u a /3 = 0) initial data. Once a (nearly) stationary solution was reached, we read off 
the value (yx) of the ut r component in the near- horizon limit (r* -C 0). Next, we evolved the 
vacuum equations again, this time with monopolar Gaussian initial data with amplitude c (and 
some width), and read off a different value (2/2) of ut r in the same limit. Finally, exploiting the 
linearity of the equations, we "tuned" the amplitude of the Gaussian initial data to —cy\j (j/2 — Vi) 
(keeping the Gaussian width fixed), to obtain a stationary solution in which ut r = at the horizon. 
Rather than performing a third evolution, we can instead make an appropriate linear combination 
of the first and second stationary solutions. It turns out that the resulting solution is in fact static 
(up to numerical error), in the sense that ut r = everywhere and not just at the horizon. To see 
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FIG. 14: Stable GLG evolution of the m = mode. These plots show numerical data (at a suitably late 
time t max = 250M) for the m = mode from the "stabilized" evolution scheme using a generalized Lorenz 
gauge and analytic initial conditions for I — and I = 1 modes (see text). The left plot shows the radial 
profile in the equatorial plane, and the right plot shows the angular profile at r = r$ = 7M. In both plots, 
the central worldtube is visible as a trough, inside which we show u^g (rather than u£p). Note that, since 
u tr vanishes with increasing resolution (see Fig. |15| below), the data is both static and in Lorenz gauge (up 
to numerical error), as desired. The numerical solution is globally- regular but has the wrong mass-energy; 
this is easily rectified by adding a specific analytically-given monopole solution as described in the text. 
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FIG. 15: Numerical data for the tr component of the m = MP, for various grid resolutions Ar* = M/n. 
Referring to the evolution shown in Fig. [l4j the plot shows the residual non-staticity left over in the m = 
mode. It illustrates that, as the grid resolution is improved, the residual non-staticity and the deviation 
from Lorenz gauge both vanish: H a oc Ut r — > as n — > oo. 

why, consider the t component of the monopole part of the gauge constraint, which reads 

u' tr - u tt + f(r)utr = 0, (138) 

where f(r) is a certain function whose specific form is unimportant here, except the feature that its 
value is bounded at the horizon. In the stationary limit we have iitt = 0, and it is easily seen that, 



with the boundary condition ut r = at the horizon, the only solution to Eq. ( 138 ) is the trivial one, 



utr = 0. In Sec. VA4 we argued that the static, globally-regular Lorenz-gauge monopole solution 



is unique; hence it seems we have a numerical approach which recovers the desired solution (up to 
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a trivial mass correction). It remains to be shown that this method can be applied to the Kerr 
problem. 



B. The m = 1 perturbation 

Despite much experimentation with choices of H a , we have not yet found a GLG that can lead 
to a stable evolution in the I = 1, m = 1 even-parity sector. Finding such a gauge remains a high 
priority for future work. 

Recall the m = 1 numerical instability takes the form of a linearly growing (Lorenz-)gauge 
mode. For circular orbits, we can remedy the situation in an ad-hoc manner, simply by applying 
a "frequency filter" 

(m=l) 1 (m=l) fl3g s 

to our numerical solutions. This eliminates any linearly growing gauge modes, while retaining the 
oc exp(±iQt) part of the dipole perturbation, which has the physically desired time dependence. 
The main drawback of this procedure is that we must compute the derivatives numerically, and this 
has the potential to amplify numerical error. In practice, we find that the trick works reasonably 
well for orbits in the strong field, but decreases in effectiveness at larger radii, as the orbital 



frequency decreases. Figure 16 shows the effect of applying this filter in the case tq = 7M. An 
advantage of this technique is that it can also be applied to circular orbits in Kerr spacetime. 
However, it should be viewed as a stopgap until one can find a stable gauge for m = 1 evolutions. 

The limitations of the frequency-filter method are revealed when we attempt to compute the 
dissipative component of the GSF. The mode i^ m_1 ) arises primarily from an I = 3 contribution, 
since I = 1 is nonradiative. However, we find that an I = 1 contribution does arise generically 
in finite-resolution simulations in 2+1D (its amplitude reduces to zero only in the limit of infinite 
resolution). Unfortunately, for realistic resolutions we find that this spurious contribution in fact 
dominates the numerical value of . 



The problem is illustrated in Fig. 17 The figure shows filtered data for the m = 1 temporal 
and radial components of the GSF, as a function of simulation time t, for various grid resolutions 
Ar* = M/n, A9 = 7r/(8n). The plot for i^ m_1 ) shows the imprint of residual oscillations, at 
frequency £1. At low resolutions (e.g., n = 4), it is clear that the oscillations begin to grow at 
late times. At higher resolutions, the amplitude of the oscillations is reduced. Nevertheless, at 
achievable resolutions, the oscillations are still large in comparison with the magnitude of F£ , 
which makes the extrapolation to n — > oo rather fraught (even after a suitable averaging procedure), 
and leads to substantial loss of accuracy. The plot for F r shows similar oscillations, with a similar 
absolute amplitude. In this case, however, the oscillations are very small in comparison with the 
magnitude of i^" 1-1 ^. 

The above limitations emphasize the need for a more systematic stabilization method for the 
m = 1 mode. We will revisit the problem in Paper IV of the series, in the Kerr context. 



VII. NUMERICAL RESULTS: PART II 

Let us now present a selection of numerical results for the total GSF, and compare with previous 
results available in the literature [2T1 EZ] • 

The accuracy of any numerical calculation is naturally limited by various sources of error. In 
Sec. IV of Paper I, and Sec. IVF of Paper II, we described three important sources: discretization 
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t/M 



FIG. 16: Application of a frequency filter, Eq. (1391, to the m = 1 mode. Here we show numerical data 
on the worldline for tq = 7M, for (a) the unfiltered field components u^p, whose envelops exhibit a linear- 
in-i growth; (b) the components after filtering, which appear stationary; and (c) the radial self-force mode 
Fr m ~ , after the filter has been applied. The m = 1 mode contribution overwhelmingly dominates the total 
radial GSF, whose value is - 0.0301(^/M) 2 (cf. Table "nl. 



error, due to the use of a grid of finite resolution; relaxation error, due to the residual effect of 
imperfect initial data; and m-mode summation error, due to approximating an infinite sum with 
a finite number of terms. In principle, these errors can be reduced by increasing the resolution 
(decreasing Ar*, AO); running for longer (increasing t max ); and computing more m- modes (increas- 
ing m max ), respectively. In practice, the situation is more subtle, since the sources of error are 
co-dependent. For example, computing more modes requires additional resolution to resolve the 
sharper physical features and faster oscillation (~ exp(— imVtt)) of the high-m modes. We applied 
the methods of Papers I & II to reduce these errors as far as possible (i.e. Richardson extrapolation; 
power-law fitting for late-time tails; and fitting the large- m. to an analytic model). The residual 
errors were then estimated in the standard way, by comparing results from different resolutions 
n max , times i max and high-m cut-offs m max . 

To be specific, to obtain the total GSF for a given orbital radius ro, we computed all modes 
up to m max = 15 for resolutions n = 4,6,8,12 and 16, with Ar* = M/n and A9 = 7r/(8n). 
Since the m- modes ^flif fall off* rather slowly, as m~ 2 , at large m, it is important to estimate 
the contribution from the remaining m > m max modes by fitting to an analytic model, ^Ijf ~ 

A/m 2 + B/m 3 + C/m 4 . On the other hand, the m-modes fall off exponentially fast with 

m, so for this component the large-m contribution is negligible — cf. Figure [9j For the modes 
m = and m = 1 we used the method described in the preceding section. As discussed above, 
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FIG. 17: Filtered numerical data for Fj} 71 1 '. The plots show the temporal and radial components of the 
GSF, j^™ -1 ) and i<^ m_1 \ extracted on the worldline as a function of time, for r — 6M and with various 
resolutions Ar* = M/n, A9 — 7r/(8n), after the application of the frequency filter (139). Both plots exhibit 
non-decaying oscillations that persist after filtering (which removes the linear- in-i instability). In the case 
of F t (left panel) the oscillations are large in comparison with the accurate value, which is quoted here from 
Ref. [104] and shown as a dash-dot black line. In the case of F r (right panel) the oscillations lead to a much 
smaller relative error. 



the m = 1 contribution to suffers from a large numerical error, especially at large rp. To 
obtain a final value for F^ eU (r > 10M), we allowed ourselves here to use the crude estimate 
p(m—i) ^ 25Af p( m - 2 ) f rom j^ e f |iQ3j - The error in m = 1 (due to the limitation of the frequency- 
filter method for eliminating gauge modes) is more perfidious than a relaxation error, because it 
does not necessarily diminish as t max increases (see Fig. 17). We found that it dominates the error 
budget for our calculation of F\ f . However, for calculations of and H the discretization 
and m-mode summation errors are dominant, and broadly comparable in magnitude at these 
resolutions. 

Table [TT] displays an indicative sample of results for the total GSF, as well as results for the 
Detweiler MP invariant H defined in Eq. (22). The error estimate, shown as a digit in parantheses, 
is found by combining the various error estimates in quadrature. Comparison with the accurate 
results of Refs. |21t 157] (obtained via a 1+1D implementation) shows that our method recovers the 
GSF and H to within a few parts in 10 . In all cases, our results agree with those of [211 [57] to 
within the estimated error in our data. 



VIII. SUMMARY AND OUTLOOK 

This paper — the third in our series — represents a significant step towards the goal of computing 
the GSF on the Kerr spacetime. Let us take this opportunity to review the achievements herein, 
and the challenges that remain. 

In this work we have reported on the first implementation of the m-mode scheme for a practical 
GSF computation. The implementation represents the first realisation of the approach outlined in 
2005 in Sec. V of Ref. [T7]. In addition to demonstrating that the GSF may be computed with a 
2+ ID time-domain implementation utilizing the puncture formulation, we also have shown that the 
GSF may be computed entirely in the time domain (at least for circular orbits). That is, whereas 
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Time-domain GSF results 






(M 2 /^ 2 )F t sclf 


(M 2 /n 2 )F^ li 


r = 6M 


-5.2355(6) a 
-5.23602 


1.3299(2) 3 
1.32984 


3.6695(7) 2 
3.66992 


r = 7M 


-4.0314(4) a 
-4.03177 


5.293(2) 4 
5.29358 


3.0096(3) 2 
3.00985 


r = 8M 


-3.3022(3) a 
-3.30239 


2-482(2) 
2.48055 


2-4475(4) 
2.44769 


r a = 10M 


-2.4462(6) a 
-2.44630 XiU 


7.348(7)* 5 
7.35254 XiU 


1.6736(2) 2 
1.67369 


r = 14M 


-1.6270(4) _ 
-1.62705 XiU 


1.2583(9)* 5 
1.25872 XiU 


9.0685(3) 3 
9.06858 


r = 20M 


-1.0889(3) ! 
-1.08893 


2.033(4)* e 
2.02994 


4.620(1) 3 
4.61896 XiU 



TABLE II: Numerical results for the temporal and radial components of the GSF in Lorenz gauge, F t self 
and F, sclf , and for the Detweiler MP invariant H [defined in Eq. (22 1], for a range of orbital radii fo- 
The table compares results from our 2+ ID m-mode implementation (upper entries), with results from the 
1+1D Z-mode implementation of Barack and Sago [21] and Barack, Sago and Detweiler [57] (lower entries). 
Parenthetical figures in the upper values indicate the estimated error bar on the last quoted decimals; in 
the lower entries all figures are significant. The values of marked with an asterisk were computed 
by replacing the m = 1 contribution (which suffers from large numerical error) with the crude estimate 



(m=l) 



25AfF t (m_2) /(896r ) from Ref. 



In the large-ro limit, H, F* cii and F 4 sel£ scale as 



and r 5 , respectively, which explains the relative loss of accuracy in the dissipative component F t . All 
values were computed from runs with resolutions Ar* = M/n, A9 — n/(8n) where n = 4,6,8, 12, 16, and 
W = 300M. 



previous studies resorted to frequency-domain methods to obtain the low multipoles, / = 1 and 
I = 0, we have persevered to tackle the problem of obtaining these modes within a time-domain 
scheme. 

Our approach is based on a constraint-damped Z4 formulation of the linearized Einstein equa- 
tions, with the constraints being provided by the Lorenz-gauge conditions. In Ref. [T7] it was shown 
that this approach works well within the 1+1 D setting, for modes I > 2. Here, we have shown for 
the first time that the approach also works in the 2+ ID setting, for modes m > 2, provided that we 
make a judicious choice of constraint damping. Furthermore, in our formulation, the modes m = 
and m = 1 are found to be free from exponential-in-t instabilities. Nevertheless, these modes are 
still plagued by linear-in-t instabilities, which are associated with the low-multipole I < 2 sector. 
In Sec. [V] we argued that the linear- in-t growth is attributed to certain (Lorenz-) gauge vacuum 
modes, which are globally regular and have the correct behavior on the future horizon. Such gauge 
modes are clearly "unphysical" (as they are inconsistent with the helical symmetry of the perturbed 
spacetime) and must be eliminated from our simulations. Unfortunately, this cannot be achieved 
via application of boundary conditions alone, because (at any finite time) the undesirable gauge 
modes satisfy the same physical boundary conditions as the physical solution. Instead, we have 
been forced to reconsider the foundations of the Lorenz-gauge formulation itself. 

In Sec. VI we investigated the idea of employing a Generalized Lorenz Gauge to promote 
stability, and we identified a class of gauge drivers which could render the axisymmetric (m = 0) 
mode stable, in the sense that arbitrary initial data evolve towards a stationary {dth a p = 0) 
solution. We picked a specific GLG with the special property that a static MP (dth a p = and 
hu = 0) would automatically be in Lorenz gauge. Then we demonstrated that (in some sense) there 
was only one degree of freedom remaining in the class of stationary solutions, so that it became 
straightforward to find initial data which led us to the static, Lorenz gauge solution. 
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Thus far we have not found it possible to identify a GLG which promotes stability in the 
even-parity I = 1 sector. Instead, we have been forced to employ an ad-hoc "trick", whereby a 
frequency-filter is applied to the Lorenz-gauge evolution to eliminate the undesirable time-growing 
gauge modes. We demonstrated that this trick works sufficiently well, insomuch as it allows us to 
compute F r and H to an accuracy of a few parts in 10 4 , for strong- field circular orbits. Nevertheless, 
our handling of the m = 1 mode remains the weakest point of our analysis, and the largest source 
of numerical error. 

Let us now contemplate the obstacles that remain to impede our pursuit of the first GSF calcula- 
tion in Kerr geometry. On Kerr spacetime, we are fundamentally limited by the lack of separability 
of the tensorial equations; indeed, this was the primary motivation for the development of the m- 
mode approach. Inevitably, non-separability means that we have a more limited understanding of 
the low-multipole sector. For example, since we cannot separate the MP into I modes, we cannot 
obtain an analytic solution for the Lorenz-gauge "monopole" piece (or even define what that piece 
is). However, the analytical techniques developed in this work can take us a long way toward 
being able to construct physical mass and angular momentum perturbations in Kerr, in the Lorenz 
gauge. Building on a key theme of Sec. |Vj we may seek a basis of homogeneous Lorenz-gauge 
solutions, by considering gauge vectors satisfying = 0, as well as mass perturbations akin to 
the "conformal" solution of Eq. ( 99 ) . The energy and angular momentum content of such solutions 
may be assessed by the method of Sec. HE Armed with a sufficiently large basis of homogeneous 
solutions, we may apply post-hoc corrections to our numerical results, in order to seek a globally 
regular solution which has the correct physical content (i.e. energy and angular momentum). Once 
more, an insistence upon Lorenz gauge will lead us to a MP which is not asymptotically-regular in 
the tt component (and possibly in t(j>). 

Inevitably, we are left with some questions that can only be answered by attempting a Kerr 
implementation. Are linear-in-t gauge modes present on Kerr? We presume yes, since the scalar- 
type gauge modes (which are straightforward to find on Kerr) are implicated. Does the GLG 
approach, generalized to Kerr, restore stability to the m = mode? If so, is the class of stationary 
solutions just as simple, in the sense described above? 

Let us now outline some areas for future research. The lack of stability in the low modes strikes 
to the core of the feasibility of the time-domain Lorenz-gauge approach. Addressing stability is 
a high priority, not least because (we believe) the issue will also affect the 3+ ID time-domain 
approach of Vega et al. In 3+1D one is not so easily able to isolate the m = and m = 1 modes 
from the rest of the system, and the problem may well hinder any naive attempt to evolve the 3+1D 
system. We may seek lessons from other approaches to the Z4 formulation in Numerical Relativity. 
Using a range of (usually implicit) gauge drivers, various authors have shown that stable evolutions 
may be achieved, even for the non-linear problem. One key difference is that their approach is 
usually based around a Cartesian grid and the harmonic gauge, rather than spherical coordinates 
and the Lorenz gauge. A further complication is that moving to more general gauges will require 
some careful thought about how regularization of the GSF may be achieved and how physically- 
meaningful quantities may be extracted |61| . Nevertheless, there is clearly much scope for a fertile 
exchange of ideas. 

There are a number of ways in which the efficiency of our numerical scheme could be improved. 
First, the method of hyperboloidal slicing, in combination with compactification of the domain 
[101J, shows much potential to dramatically reduce the computational burden for time-domain 
calculations. Second, the method of mesh refinement under development by Thornburg |26|, [27] is 
ripe for application to this problem. Third, the fourth-order puncture formulation of Wardell et 
al. [33] may be easily adapted to our m-mode scheme. Fourth, a careful optimization of the initial 
data could help us minimize the seeding of junk radiation that decays away only slowly in the low 
multipoles. 
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The m-mode GSF program has four immediate priorities: (i) to find a GLG which renders the 
m = 1 mode stable; (ii) to find a GLG which is appropriate for generic orbits (bearing in mind that 
we want the transformation from the GLG to Lorenz gauge to be regular); (iii) to implement a 
fourth-order puncture scheme, which will yield both a more rapid convergence with grid resolution 
and a more rapid convergence of the m-mode sum, leading to higher numerical accuracy; and (iv) 
to press ahead with the an implementation in Kerr. In the longer term we have two aims: (i) to 
develop the method up to the point where it can handle self-consistent orbital evolutions, in which 
the influence of the GSF is applied at every time step to correct the orbital trajectory, and (ii) to 
incorporate GSF effects which are second-order in fi/M, based on the recent formulation by Pound 

PE9E5I and Gralla [53 • 
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Appendix A: Key quantities in the 2+1D field equations 



We give here explicit expressions for the matrices and -M-^g appearing in Eq. (48) and 



(51), respectively. Recall u denotes differentiation with respect to t, and vf denotes differentiation 
with respect to the tortoise coordinate r* defined in Eq. ( 39 ) . For brevity we set M = 1 throughout 



this appendix; the correct factors of M can be easily recovered, if needed, by considering the 



dimensionality of the various terms. Also for brevity we omit the suffix (m) in u 



The coupling terms in the original 2+1D field equations (48), prior to the imposition of gauge 
constraint damping, are given by 



(m) 



j^im) = 2 {2r\utr 



+ U tt - Urr) 4/ (u tt -U rr ) 2f 2 (u 9e + U 4 



+ 



(Al) 



~ ( m ) _ 2/ 2 (cos Bute + imutcf.) 2(u u + u rr - 2u' tr ) 2f 2 (u tr + d e u te ) . . 

JVl t — - | ^ a ' \ > 

A^M _ f( Uts + 2imcos 0u t<t>) , 2(n r g - u' te ) f[(4 + r)u t e + 2rd e u tr \ _ f 2 u t e ... 
M te ~ r 2 s[n 2 Q + r 2 + r 3 r 2 ' ^> 

.-Am) _ f(u t4 , -2im cos 6u tg ) 2fimu tr 2 ( u r<t> ~ u' t(j} ) f(4 + r)u t<t> f 2 u t(f> 
v r z sin 6 r A sin a r z r 6 r z 

^(m) 4f 2 (cos9u r9 + imu r ^) 2[2r 2 (u tr - u' rr ) + u„ - u tt ) 4f(u tt -u r r) 
2f 2 (2ru rr + u e e + u^ + 2rdeu re ) | 2f 3 (u e9 + u^) 

^ 1 ^2 ' y Ai> > 
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-Am) _ f(u r + 2im cos 9u r(t> ) 2f 2 [cos6(u 99 - u^) + imug^] 2(u te - u' r9 ) 



M K L = - - H ^ 

re r 2 sin 2 6* r 2 sin r 2 

| /[(4 + r)u rg + 2r9 e u rr .] _ f 2 (5u rd + 2d e ugg) , Ag . 

i^*3 y>2 



^(m) _ /(u r - 2zm cos 6>M r g) 2/[2/ cos Oug^ + im(/^ - u rr )] ^ 2(tt^ - ^) 
| /(4 + r)n r0 f 2 (5u r(f> + 2dgu e( ],) 

a^M _ 2/[^gg - + 2zmcos(9tte0] 2(n ff - u rr ) 2f(u rr + ng e + 2d e u rS ) _ 2f 2 (u ee + 



r 2 sin 2 6 1 r 3 



r 



2 



(A8) 



-( m ) _ 2/[2-ug^ - imcos8(u ee - u^)} 2/(cos 9u r ^ - imu r e) | 2/^^ + d e u ri p) 

^9<b ~ 9-2/1 2 ■ a ' 2' 

v r 2 sin 6> r z sin r z 

A -,(m) _ 2/(tt ee - + 2im cos 9u e<j ,) 4/(cos 6>-u r g + imu T(j} ) 2(u tt - u rr ) 2f(u rr + u^) 



r 2 sin 2 # r 2 sin# r 3 r 2 



2f 2 (u ee + u 4 



The components of the m-decomposed divergence = /i^J 1 ' are given by 



(A10) 



~M - -Uv f (cos 6u te + imu t(f> ) u tt - u' tr f(u tr + d e u te ) 

= A* J z t — o • a 1 2 ' l Aii J 

r z sm u r r z 

grim)- -1*2-7 _ f(cOs8u r g + imn^) U tr - | Mff - U rr ( /(tt rr + 0gU r e) ^(uge + U^) 



y2 g2j]_ $ ■y* y*3 ^>2 y>2 

(A12) 

7 (m)_ -if f/ \y _ f[cos9(u ee - Ufa,) + imue^] u t e - u' rd f (2u r8 + dgu ee ) , A1 „x 

= A* U/ r ) Zl 6— <T^ 7) I o ' V Ai,j i 

r sm 6/ r r z 



4 m) - /^[//(rsinfl)]^ = /(acosg^ + imu^) _ ug - + /(2^+jg^) 

The factors of / are introduced so that the components Z a are regular at the event horizon — see 
the discussion in the main text. 

The coupling terms in the eventual 2+1D field equations (51), after the imposition of gauge 
constraint damping, are given by 

j^(m) _ 4f (cos 6u te + imuty) 2 [u u - u rr ] 4 [u tr + u' tr - (u tt + u' tt )] 
tt r 3 sin0 r 4 r 2 

| <jf(u tt + U tr ~ U rr + 9gU W ) | 2f 2 (u dd + U^) 



M 



( m ) _ 2f [cos 0((r - 3)u t e - u r8 ) + im((r - 3)u t(j) - u ri p)] 2(u a - u„) 2[u rr + u' rr - (u tr + u' tr )] 



tr 



r 3 sin ( 



+ 



+ 



2f(u rr + u tr + dgu r 9 + d 9 u t e) 2f 2 (ru tr + rdeuffi + u ee + u. 



+ 



(A16) 
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j^(m) _ fjute + 2im cos 6u t(j) ) | 2/[cos 9(u ee - u H ) + imu e<f \ | 2[u rd + u' rB - (u t e + u' tg )} 
te r 2 sin 2 6 r 3 sin r 2 
t /[(4 + r)tt f g + Au t8 + 2a e ^g + 2rd e u tr ] f 2 u te 

.~,(m) _ f{u t ^-2im cos 6u t e) 2 f [2 cos Ouq$ + im(ru tT + u^)} 2 for0 + ^ - + u' t(f) )} 



/((4 + r)u t( p + 4n r + 2d e ue<p) _ f 2 u t $ 

^3 ^2 



(A18) 



~ (m) _ 4/(r - 3)(cos 6u r g + imu rf p) 2{u u - u rr ) kf[u u - 2u„ - d d u r e] 
r" 5 sm 6* r 4 
2/ 2 [3(uee + u^) + 2r(u rr + ^-u^)] 2f 3 (u ee + u^) 



y 3 y«2 



(A19) 



^(m) _ /(Mr6» + 2im cos 9u r(j> ) 2/(r - 3) [cos 8(u ee - u^) + imu e 



r6 r 2 sin 2 8 r 3 sin ( 



/((8 + r)n r . e + 2^^ + 2rd e u rr ) f 2 (5u r9 + 2d e u ee ) 

H 5 R j (AZUJ 



~( m ) f (u ri p — 2im cos 9u r $) 2/(2(r — 3) cos flit^ + im(r — 3)tt^ — imru rr ) 
r< t> = r 2 sin 2 r 3 sin0 



| /((8 + r)u r4> + 2d^) / 2 (5u T . (/ , + - ^, 



Appendix B: Elliptic integrals in the puncture and effective source 

In this appendix we list relevant elliptic integrals that go into the construction of the m-mode 



puncture function [Eq. (69)] and effective source [Eqs. (75) and (76)]. Below, elK(-) and elE(-) 



denote complete elliptic integrals of the first and second kinds, respectively, defined by 

fn/2 

elK(fc) = / (l - A; 2 sin 2 xT 1/2 cfa, (Bl) 
Jo 

e\E{k)= (1 - k 2 sin 2 x) 1/2 dx. (B2) 

Jo 
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/r = 



jm 
1 2 



The relevant integrals are 

C = £ e~ m54, d(H) = [pffif (p)elK(7) + P^(p)dE( 7 )] , (B3) 

£ t~T e- m54, d(H) = ^ te(p)elK( 7 ) + p- 2 ^(p)elE( 7 )] , (B4) 

£ e v 3 cos 5<Pe-^d(5<P) = ^ fc(p)elK( 7 ) + p- 2 p^(p)elE( 7 )] , (B5) 
/ 3 m £ e^ 5 cos 2 (50/2) e-^d(^) = fe(p)elK( 7 ) + ^ 2 p^(p)elE( 7 )] , (B6) 

£ e v 5 sin 2 (50) e- imS *d{8<t>) = ^ [ P S,(p)elK( 7 ) + p- 2 p^(p)elE( 7 )] , (B7) 

£ e p 5 sin 2 (50/2) e-^dm = ^ fe(p)elK( 7 ) + /T^/Wt)] , (B8) 

C = £ ^ 3 sin 2 (50) e- imS *d{8<t>) = ^ »(p)eIK( 7 ) + P^(p)elE( 7 )] , (B9) 

/ 7 m = £ e^ 1 ccnS* e- im5 <?d(6<l>) = ^ fe(p)eIK( 7 ) + p^(p)elE( 7 )] , (BIO) 



Tin 
1 A 



Tin 
2 5 



Here 



with 



and 

-1 ■„ r jl -im8<b j/rj.\ ~ l P 



JT = J e p 1 sin50e- jm5 M^) = ^teelK(Vp) + C^lE(Vp)], (Bll) 

JT = £ ep 3 sin50e-^d(50) = [^elK(i/p) + p 2 q? E e\E(i/ p)] , (B12) 

J 2 m = £ e p 3 sin 50 cos 50 e - m<5 <^(50) = y^J fe"kelK( 7 ) + ^elE( 7 )] , (B13) 

JT = £ ^ 5 sin 50 cos 2 (50/2) e~^(50) = ^ [g&elK( 7 ) + p^elE^)] , (B14) 

JT = £ ^ 5 sin 50 sin 2 (50) e - imS *d(6<l>) = [^elK^/p) + p 2 ^elE(i/p)] , (B15) 
JT = £ ^ 5 sin 50 sin 2 (50/2) e - m<5 ^(50) = ^ feelK(i/p) + p 2 g 5 >lE(Vp)] . (B16) 

p 2 ^ A/(4B), 7 ^(l + p 2 )- 1/2 , (B17) 

^ = P rr 5r 2 + P e0 5# 2 + Q rr 5r 3 + Q 0e 5r5# 2 , B = P^ + Q H 5r, (B18) 



PrT = /o _1 ! ^=r 2 , = r ° 3/ ° , (B19) 

ro — oM 



Qrr = -^, Qee = r , Q44 = r ( j^- ^ ) . (B20) 
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™ Pgk(p) 



1 | (4p 2 + l) (4p 2 + 3) 

2 | (2p 2 + 1) (64p 4 + 64p 2 + 5) 



q 8192 n 8 _|_ 16384 ^6 , 6176 „4 , 6304 „2 , n 

J 35 P 35 P "T" 21 P + 105 P + Z 

4 gfg (2p 2 + 1) (81920p 8 + 163840p 6 + 97408/ + 15488p 2 + 315) 

r 1048576 12 , 1048576 10 , 17686528 8 , 9158656 fi , 427744 4 , 180832 9 
231 l 77 g 1155 P '■' 1155 g 231 P 1155 g 



m p™ E {p) 



1 -f (2p 2 + l)(p 2 + l) 

2 -§ (p 2 + l) (64p 4 + 64p 2 + 9) 

3 (2p 2 + 1) (p 2 + l) (768p 4 + 768p 2 + 53) 

4 -gfg (p 2 + 1) (81920p 8 + 163840p 6 + 102528p 4 + 20608p 2 + 683) 
5-rlfg (2p 2 + l) (p 2 + l) (l63840p 8 + 327680p 6 + 194304p 4 + 30464p 2 + 659) 



m p? K (p) 

4p 2 + 2 

1 fp 4 + fp 2 + 2 

2 A (8p 2 + 5) (2p 2 + 1) (8p 2 + 3) 

o 2048 n S , 4096 n 6 , 1712 d , 2416 _2 , 9 
35 P ~r 35 P + 21 P 105 1° + Z 

4 gfg (2p 2 + l) (l6384p 8 + 32768p 6 + 21632p 4 + 5248p 2 + 315) 

r 524288 12 , 524288 10 , 9181184 8 , 5255168 6 i 100624 4 , 192496 2 

_ 693 P 231 g 3465 P r 3465 ^ : 231 g : 3465 ^ 



m 


1 

2 
3 

4 
5 



Pt"e(p) 



(2p 2 



lg5 
315 

'3465 



2p 

P 2 + 
(2p 2 



1) (P 2 + 
1) (64p 4 

+iW 

- 1) (I6384p 8 
1) (P 2 + 1 



64p 2 + 19) 

1) (384p 4 + 384p 2 + 79) 

- 32768p 6 + 22656p 4 + 6272p 2 + 523) 
(l63840p 8 + 327680p 6 + 215424p 4 + 51584p 2 



3079) 



TABLE III: The polynomials p™ K and p™ E appearing in Eqs. (B9) and (BIO), for n = 6, 7 and m = 0, .. . ,5 



We remind /o = 1 — 2M/rQ. The quantities p^- and in the above expressions are polynomials 
in p 2 . Some of these polynomials (for m = 0, . . . , 5) were tabulated in Tables I and II of Ref. |81j 
(where the notation was used in place of our P^/g), and a few others in Table VIII of Paper 
I. All remaining polynomials relevant for the current work are given in Tables |Hlj |IV| and |V| below, 
again for m = 0, . . . , 5. Higher-m polynomials can be similarly obtained using a symbolic algebra 
package. 



Appendix C: m-mode decomposition of the effective source 



In Sec. IIIC3 we introduced the effective source, and in Eq. (74) we gave an expression for 



its m-mode decomposition. Here we complete the description by giving explicit expressions for 
the coefficients and d k a ^ featuring in Eq. (76) for AZ^ . Here O is the orbital frequency, A r , 
Aq, B r , Bg, etc., denote the partial derivatives of A and B defined in Eq. ( B18| ), and we define 
E a p = u a up + D a p5r where D a p were given in (62)-(66). After setting M = 1 the (m- independent) 
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m Qokip) 



1 -f (/ + l) 

2 -lM(2/ + l)(/ + l) 

3 -|f (/ + l) (128/ + 128/ + 27) 

4 -f^(2/ + l)(/ + l)(32/ + 32/ + 5) 

5 -||(/ + l) (32768/ + 65536/ + 44160/ + 11392/ + 875) 





1 §(2/ + l) 

2 § (l6p 4 + 16/ + 1) 
2/ + 1) (128/ + 128/ + 3) 

i o5536 n S I 131072 n 6 , 26624 4 , 2048 _2 , 32 

^ 315 ^ " r 315 P " r 105 ^ ^ 45 ^ " r 63 

5 (V + y (32768/ + 65536/ + 38016/ + 5248/ + 35) 





1 2(2/ + l) 

2 | (V + 1) (4/ + 3) 

3 § (2/ + 1) (128/ + 128/ + 15) 

a 16384 n 8 , 32768 fi , 12800 4 , 14848 Jl , o 

4 Tag-** + ^5~^ + T05^ + 8 

5 ^ (2/ + 1) (32768p 8 + 65536/ + 40576/ + 7808/ + 315) 



m iTe(p) 



1 -4 

2 -f (2/ + 1) 

3 _512p4 _ 512 ^2 _ 92 

4 -f|g (2/ + 1) (96/ + 96/ + 11) 

r. _ 131072 8 _ 262144 6 _ 18944 4 _ 5632 2 _ 2252 
63 P 63 P 7 P 9 P 63 



m q™ K (p) 



1 § (4/ + 1) (4/ + 3) 



2 § (2/ + l) (32/ + 32/ + 5) 

o 8192 s , 16384 ^6 i 2208 _4 , 2848 _2 , o 

35^ i "35/ 5i "7/ ,i "35/ 5+D 

4 Jg (2/ + 1) (20480/ + 40960/ + 26368/ + 5888/ + 315) 

r 1048576 _12 , 1048576 „10 j_ 3645440 s , 2048000 „6 568544 „4 , 22944 „2 , in 

231 P 77 P "' 231 P " 231 231 P ~~ ^77~^ + 1U 



m 



1 

2 
3 
4 
5 



q?E(p) 
o 

-¥(2/ 



16 

231 



2/r 



+ 1) (/ + !) 
1) (32/ + 32p 2 + 7) 
+ 1) (/ + 1) (256/ + 256/ + 41) 
f 1) (20480/ + 40960/ + 27648/ 
+ 1) (/ + 1) (32768/ + 65536/ 4 



- 7168/ - 
42240/ 



533) 
■ 9472/ 



519) 



TABLE IV: The polynomials q" n l K and q™ E appearing in Eqs. (Bll ) - (B13), for n = 0, 1, 2 and 
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m q? K (p) 



1 + 

2 -|(32p 4 + 40 i o 2 + ll) 

3 -i (8/9 2 + 5) (256p 4 + 30 V + 63) 

4 _ 8192 8 _ 96256 6 _ 77824 4 _ 8256 2 _ 70 

21 " 105 " 105 P 35 P 3 

5 (8p 2 + 1) (32768p 8 + 90112p 6 + 89344p 4 + 37520p 2 + 5565) 



m q™ E (p) 



1 i(8 P 2 + l)(p 2 + l) 

2 Kp 2 + 1) (32p 4 + 24p 2 + 1) 

3 hip 2 + l) (2048P 6 + 2688p 4 + 808P 2 + 15) 

4 jL(p 2 + 1) (20480p 8 + 37888P 6 + 21248p 4 + 3488p 2 + 35) 

5 (|Q 2 + 1) (262144p 10 + 622592^ + 509952/q 6 + 165248p 4 + 17080p 2 + 105) 



m q? K {p) 



1 I (4p 2 + 1) (4p 2 + 3) 

2 j| (2p 2 + 1) (I28p 4 + 128p 2 + 15) 

o 8192 ^8 i 16384 „6 , 10592 „4 i 800 „2 , c 
° 21 ^ "r" 21 P "T" 21 ^ "r" 7 P +0 

4 gfg (2p 2 + 1) (40960p 8 + 81920p 6 + 50048p 4 + 9088p 2 + 315) 

r 1048576 12 , 1048576 10 , 24977408 8 , 13254656 6 , 1091936 4 , 318496 2 
_ 99 P 33 P 1 693 " 693 g 231 " : 693 " 



m qf E (p) 



1 -f (2p 2 + l) 

n _ 1024 n i _ 1024 „2 184 
15 ' 15 ' 15 

3 -if(2p 2 + l) (256p 4 + 256p 2 + 27) 

4 _ 131072 8 _ 262144 6 _ 280576 4 26624 2 _ 9328 

63, ^ 63 " 105 " 45 " 315, 

(2p 2 + 1) (229376p 8 + 458752p 6 + 278784^ + 49408p 2 + 1697) 



™> qsk(p) 



1 i(p 2 + l)(8p 2 + 3) 

2 | (p 2 + 1) (32p 4 + 24p 2 + 3) 

3 % (p 2 + 1) (8p 2 + 3) (256p 4 + 208p 2 + 15) 

4 (p 2 + 1) (20480/9 8 + 33792^ + 17408p 4 + 2976/0 2 + 105) 

5 X (p 2 + 1) (8p 2 + 7) (32768p 8 + 40960p 6 + 15616P 4 + 1904p 2 + 45) 



m q™ E {p) 



1 -I(8p 2 + 7) 

2 -|, 4 -f p2 - 3 

1024 fi 576 A 788 „2 51 



13 15 P 5 P 15 P 10 

8704 ^4 _ 

r .'; t072 „1Q " Igojj „8 15 320512 „6"~" 125248 "3 _ 9 Q o„2 _ 1231 
63 g 3 g 63 g 63 g AVL P 126 



i _ n 8 _ 88064 6 "8704 4 tl 14528 2 258 



TABLE V: The polynomials q"^- and g^g appearing in Eqs. (B14) - (B16 1, for n — 3, 4, 5 and m 
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coefficients c^g and 4/3 read 

o 2(r - 1) 2(r-3) 2/ 

% - r4/ r 2 Dtt + r s sin 2 6 ^ { ' 

4 = {2E tt /r 2 - fD tt ) (A r + 2B r ) , (C2) 

4 = -2 (2E tt /r 2 - fD tt ) B r , (C3) 

4 = 4D (r OB/r 2 , (C4) 

4 = -4A r n/r 2 , (C5) 

4 = -8(r - 1) Ar/4 (C6) 

4 = AflBEu/ifr 2 ) + AB/ (r 4 sin 2 + 2/(4- + 2B T )D tr /r 2 , (C7) 

4 = -4fB r D tr /r 2 , (C8) 

4 = 2^/(r 4 / 2 ) - 2£^/(r 5 sin 2 0), (C9) 

4 = -25 [fiAr/(r 2 /) + ^/(r 4 sin 2 9)] , (CIO) 

4 = 2flD tr /(r 2 f) + 2D r<p /(r' l S m 2 9), (Cll) 

4= [(-r 4 ft 2 + 2r 2 - 6r + 4)/(r 4 /) + (r 2 sin 2 0)~ 2 ] Ar, (C12) 

4 = (A. + 2B r )D tr /r 2 - 2(r-3)B^/(r 4 / sin 2 0), (C13) 

4 = [-2( J Br 3 ^ 2 + rfB r )/(r 3 f) + 25/(r 2 sin 2 0)] D tr , (C14) 

4 = -2 cos 9E <t><t> /(r A sin 3 0), (C15) 

4 = -2cos^^/(r 2 sin 3 0) + /AArA, (C16) 

4 = 5iVr 3 , (C17) 

4 = A e E t4> cos #/(r 2 sin 0) + ((r - l)£^/r 2 - fD ttf> ) (A r + 2A) , (C18) 

c% = -2B r ((r - l)iVr 2 " /A*) , (C19) 

4 = -25 (/Ar/r - SlD r4> /r 2 ) , (C20) 

4 = 2 (/Ar/r - fiA*/r 2 ) , (C21) 

4 = -4(r - 1) D r4> /r\ (C22) 

4 = 2BE (j>4> /(r*sm 2 9) + f(A r + 2B r )D r( /,/r 2 + 2nBE t(f> /{r 2 f), (C23) 

4 = -2/AA^/r 2 , (C24) 

4 = -2(3r-7)A t /(r 4 / 3 ) + 2(r-5)^/(r 5 /sin 2 0), (C25) 

4 = 4(r-3) J BA^/(r 4 /sin 2 0), (C26) 

4 = -4(r-3)A<4(r 4 /sin 2 #), (C27) 

c°, = 2cos#A^(r-3)/(r 4 /sin 3 fl), (C28) 

c% = 2BD r ^cos9/(r 2 sm 3 ff), (C29) 

c 7 r6 = -2D r<p cos9/(r 2 sm 3 9), (C30) 

d% = D r4> [-(r 4 n 2 - Ar 2 + 19r - 18)/(r 2 /) + 1/ sin 2 0] /r 2 , (C31) 

4 = -A4( r - 4 )(^ + 2A)+Acot0]/r 2 -25(r-3)^/(r 4 /sin 2 0), (C32) 

d% = 2D r4> [(r - 4) A - r 3 n 2 B/(fr) + 5/ sin 2 (9] /r 2 , (C33) 

4, = 2At/(r/) + 2£7^/(r 2 sin 2 6) [-/ + 1/ sin 2 9] , (C34) 

dg = 4/cos0A*/(rsin0), (C35) 

dj = A e fD r4 ,/r - 2BE H cos #/(r 2 sin 3 9), (C36) 

cj = 2E tt sin 2 0/(r/) - 4E^/r 3 - 2(r - 3) £>^/r 2 + 2£7^/(r 2 sin 2 0), (C37) 

4 = /(2£^ -rD H> )(A r + 2B r )/r + 2A E <M> cos #/(r 2 sin0), (C38) 

= -2/(2^ - rD^Br/r, 52 (C39) 

c% = -4fBD r<p /r, (C40) 



It should be noted that the appropriate effective source depends upon the choice of gauge constraint 
damping (Sec. II C); the expressions above correspond to the choice made in Eq. (50). 



Appendix D: Projection onto a tensor-harmonic basis 



We give below formulas for constructing the multipole lm-m.ode functions huL[r,{) (as defined 

by Barack and Lousto in [17] ) from our m-mode variables u^o (t,r,9). We omit the suffix (m) for 
brevity. 



h 



(10) 
Ira 



r,t) 
r,t) 



r,t) 
r,t) 



r,t) = 2tt sw.9(uu + u rr )Y* m dO, 



r,t) = Air sin 9 u tr Y{ f m d9, 



r,t) = 2ir I sin# (u u - u rr ) Y l * m d9, 

[sin u t o d e - imu t ^} Y* m d9, 
[sin 9 u r g dg — imu r ^\ Yi m d9, 



4tt 
4tt 



r,t) = 2ir sva.9(u e e + u 4>(t) )Y l * m d9, 



r,t) = 2tt 



4tt 
4tt 



sm6(uee-u^)D 2 + 2ue ( / > Di Y l * m d9, 
-imu w - sin 9u t(t> d e ] Y,* m d9, 
—imu r Q — sin#u r 0<%] Y l ^ n d9, 



r,t) = 2tt (u e9 - u^)Di - 2 sin 9u 9( f>D2 



Y? m d0. 



(Dl) 
(D2) 
(D3) 
(D4) 
(D5) 
(D6) 
(D7) 
(D8) 
(D9) 
(D10) 



Here Y[ m (9,(fr) are the usual spherical harmonics, an asterisk denotes complex conjugation, D\ 
2 (de — cot 9) and L>2 = dee — cot 9dg — (sin9)~ 2 d ( p ( j > . 
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